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FOREWORD 


Work described in this report was performed by personnel of Textron’s Bell 
Aerosystems Company for the Structures and Mechanics Division of the NASA Manned 
Spacecraft Center, Houston, Texas under contract NAS 9-3528. The work was per- 
formed in the period September 12, 1964 to March 1, 1966 under the direction of 
Mr. R.H. Gallagher, who acted as Principal Investigator. The effort was monitored 
by Mr. F.J. Stebbins of the Structures and Mechanics Division, NASA Manned Space- 
craft Center. 

The subject effort entailed contributions from many disciplines and therefore 
represents the contributions of numerous individuals other than the cited authors. Hie 
authors are particularly appreciative of the support provided this project by Messrs. 
F. Braun, and W. Luberacki. 
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ABSTRACT 



/ i 



Analyses are performed of the Apollo Aft Heat Shield for conditions of water 
impact and thermal stress. The water impact condition is first treated as a quasi- 
static problem in the determination of localized stresses of a critical nature. Two 
separate temperature profiles are examined in the solution for thermal stresses . 

Dynamic response and elastic instability analyses are performed for a limited 
number of conditions in order to confirm the influence of these effects with respect to 
the quasi-static solution. The concepts of discrete element analysis, as they pertain 
to the Apollo Aft Heat Shield, are reviewed and more advanced and satisfactory 
idealization procedures are defined. 
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I. INTRODUCTION 


In September, 1964, the NASA Manned Spacecraft Center contracted with Bell 
Aerosystems Company for a six-month study of the application of the discrete element 
structural analysis method to the prediction of thermal stresses in the Apollo command 
module aft heat shield. Shortly thereafter, in Dec., 1964, the Apollo aft heat shield sus- 
tained damage during tank-testing, occasioning efforts directed towards structural 
design changes, particularly with respect to the scalloping of the face sheets. Bell 
Aerosystems Company was asked to assist in this work, via an extension of the con- 
tract, and numerous analyses were performed to guide the specific nature of the 
scalloping to be effected. 

The latter effort, which was conducted within a very short period of time, was 
amplified in March 1965 into a more extensive examination of the heat shield with 
respect to water impact and thermal stressing conditions. This work is described in 
detail herein. Other objectives include the substantiation for this class of structure of 
the principal analysis technique — the matrix discrete element approach — through 
comparisons with experimental and classical analysis results, and the development of 
design data for the general heat shield problem. 

To meet these objectives, the subject effort was divided into five areas of 
activity each with subsidiary work items. These are the 

(1) performance of basic design analyses for water impact, treated as a quasi - 
static problem, and for thermal stress conditions. 

(2) performance of dynamic response analyses for a limited number of 
conditions. 

(3) determination of buckling pressures for water impact conditions 

(4) examination of discrete element stiffness equations for more accurate 
matrix analysis of heat shield structures . 

(5) development of a classical solution and, from this, design charts for 
anticipated heat shield design problems. 

Section II of this report describes the work performed in connection with basic 
design analyses. These analyses exclude the effects of dynamic behavior and elastic 
instability. Included are assessments of convergence characteristics, two thermal 
stress analyses, and an extensive series of analyses representing a quasi-static 
treatment of water impact conditions. 

Analyses in which the effects of dynamic response are neglected generate results 
which are expected to be conservative. Since it is desirable to obtain a quantitative 
estimate of the degree of conservatism, if any, involved in the static analysis, the 
dynamic response of the Apollo aft heat shield to water impact was determined. One 
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impact angle was chosen and histories of the loading at all affected points were 
defined. The results of this effort appear in Section m. 


Elastic instability effects are examined in Section IV. The water impact con- 
dition produces membrane stress throughout the heat shield and it is conceivable that 
these stresses could be large enough to produce a failure of the heat shield in an 
elastic instability mode. Thus, using a discrete element approach to linear instability 
analysis, computations were performed for three distributions of pressure to deter- 
mine the intensity of pressure to cause buckling. The results, which were modified by 
means of empirical factors to account for the influence of nonlinear effects on 
buckling, indicate a small but positive margin of safety. 

The questions of an appropriate discrete element idealization for the Apollo aft 
shield are examined in Section V. The significance of the heat shield as a sandwich 
structure is considered and it is found that principal modes of behavior are represented 
with sufficient accuracy by means of conventional thin plate formulations. A discrete 
element formulation for sandwich structures is advanced, however, and it is shown that 
it holds promise for accurate representation of shallow sandwich shells, provided a 
high degree of gridwork refinement is feasible. The existing discrete element formu- 
lations for shallow shell analysis are reviewed and it is found that new approached, 
based on satisfaction of interelement compatibility requirements, hold promise for 
resolution of difficulties inherent in the present formulations. 


A classical approach, to the problem of analyzing the Apollo heat shield is pre~ 
sented in Section VI. This approach covers both thermal stress and applied load 
situations and is extended to account for the overhang of the heat shield beyond the 
bolt circle. Conclusions are presented in Section VII. 


The detailed aspects of the discrete element procedures for instability and 
dynamic analysis are given in two appendixes. Another Appendix describes the Bell 
General Purpose Structural Analysis Computer Program. 
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n. BASIC DESIGN ANALYSES 


A. OBJECTIVES 

Within the design analyses, a number of classes of problems have been treated. 
These include 


(1) Analyses to determine the effects of grid refinements 

(2) Analyses in support of the design activity which took place during January 
1965 

(3) Analyses to determine critical water impact cases. 

(4) Analyses for the stresses and displacements of the shield under elevated 
temperature conditions 

The actual structure of the Aft Apollo Heat Shield is drawn in Figure E-l. This 
is the "Block 2" heat shield, with sculptured face sheets. Since the specified loadings, 
which arise from either water impact pressures or temperature, are symmetric about 
the axis of geometric symmetry, only one half the heat shield is considered in 
analyses . 


Detailed objectives of the above three types 


rvf r» wn IttO/ao f r 

v/i aiiaxjf 0V/u , i/v 


-v+V»QT' tiri+Vi fVi c\ VOOllltc 
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obtained, will be presented following a brief review of the analysis procedure and the 
material properties common to all analyses. 


Included were an extensive series of analyses performed during early phases of 
the project in conjunction with unscalloped designs, various arrangements of the 
scalloping, and with the verification of test results from experiments conducted on the 
ATR-209 heat shield specimen. 


B. DESCRIPTION OF PROCEDURE 


The analyses of the Apollo aft heat shield described in this section were per- 
formed with use of the Bell General Purpose Structural Analysis Program. This pro- 
gram is based on the discrete element approach to matrix structural analysis . A 
complete description of the theoretical basis of the program, together with an outline 
of program capabilities, is presented in Appendix A. 

Figure II-2 illustrates a typical idealization of the heat shield structure by 
means of discrete elements. In the idealization shown, five different types of element 
are employed. 

(1) quadrilateral plate in plane stress 

(2) triangular plate in plane stress. 
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(3) quadrilateral plate in bending 

(4) triangular plate in bending 

(5) axial -flexural-twist element 

The first four elements are actually representative of two physical elements of 
the structure, the quadrilateral and triangle. These elements appear in Figure 11-3. 
From a theoretical standpoint, it is necessary to divide the deformational behavior of 
these elements into plane stress and flexural behavior, respectively. 

A detailed development of the pertinent analytical relationships for these plate 
elements, in the form of the pertinent stiffness matrices and terms for thermal stress 
analysis, is given in Reference 3. It should be noted that the flexural properties for 
both elements are based on simple assumed displacement assumptions of polynomial 
form. 


In establishing the membrane and flexural rigidities of the plate element, for 
water impact analysis, the sandwich form of construction requires special considera- 
tion. Denoting the face thickness as tf, the membrane ridigity, Dm can be written as 

D , E (A t 4 ) (n-i) 

where is Poisson's ratio. Any direct stress carrying capacity of the ablator, 
bond, and core is disregarded in the water impact computations. (The ablator is 
included in analyses for thermal stress conditions.) The membrane rigidity is 
normally expressed as 

P*v - (H-2) 

(t-A- ) 

Hence, by comparison of II-l and II-2 


t M “ A tf 

For the flexural rigidity at a plate, 0 ^ , one can write 

D k = ill 

In the case of the sandwich, however, the flexural rigidity is actually 
Comparing II -4 and n-5 


(II— 3) 


<n-4) 


(H-5) 


(II -6) 
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The following is a listing of the equivalent thicknesses, and ^ , as computed 
from equations n-3 and II-5, for the various face thicknesses for sandwich portions of 
interest (h = 2.0 in.) are: 


*f 

t 

m 

•b 

0.008 

0.016 

0.577 

0.012 

0.024 

0.661 

0.020 

0.040 

0.783 

0.030 

0.060 

0.896 

0.050 

0.100 

1.063 

* 

0.025 

0.513 


♦Corrugated torus flange 


The fifth element cited previously is the axial flexural -torsional element used to 
represent the ring at the circumference of the heat shields in the analyses where the 
corrugated torus was included. The section properties used were: 


Area 

Moment of Inertia about the X' Axis 
Moment of Inertia about the Y' Axis 
Torsional Rigidity 


0.235 in. 2 
0.309 x 10 -2 in. 4 
0.650 x 10 -2 in. 4 
0.710 x 10~ 2 in. 4 
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Points 1-4 

b. Quadrilateral Plate Element 


Figure n-3. Plate Elements for Heat Shield Analysis 
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TABLE H-l 

NODE POINT COORDINATES - COARSE GRID 
(See Figure n-2 for Layout) 


Node 

Point 


1 

-72.21 

0 

160.066 

2 

-58.100 

0 

165.710 

l 

-47.982 

0 

168.919 

4 

-29.266 

0 

173.144 

5 

-10.550 

0 

175.283 

6' 

9.570 

0 

175.338 

7 

30.170 

0 

172.989 

8 

58.100 

0 

165.710 

9 

72.210 

0 

160.066 

10 

-72.032 

5.063 

160.066 

11 

-57.879 

5.063 

165.710 

12 

-49.021 

7.550 

168.450 

13 

-31.001 

12.606 

172.382 

14 

-12.981 

17.663 

174.227 

15 

6.392 

23.094 

173.957 

16 

26.226 

28.659 

171.250 

17 

46.797 

34.433 

165.710 

18 

61.121 

38.451 

160.066 

19 

-71.210 

11.973 

160.066 

20 

-56.853 

11.973 

165.710 

21 

-52.259 

14.859 

166.984 

22 

-36.407 

24.809 

169.985 

23 

-20.555 

34.759 

170.894 

24 

-3.514 

45.456 

169.578 

25 

13.935 

56.404 

165.710 

26 

29.260 

66.016 

160.066 

27 

-67.158 

26.535 

160.066 

28 

-55.050 

18.577 

165.710 

29 

-45.854 

35.680 

165.710 

30 

-57.962 

43.066 

160.066 

31 

-32.046 

48.463 

165.710 

32 

-44.153 

57.139 

160.066 

33 

-12.055 

56.836 

165.710 

34 

-24.163 

68.047 

160.066 
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MATERIAL PROPERTIES 


The cross section of the Aft Apollo Command Module Heat Shield is shown in 
Figure n-1. Since the properties of the bond are neglected in all analyses described 
in this report, the materials of interest are the ablator, the face sheets of the sand- 
wich skins, and the core of the sandwich. Only the face sheet properties play a role 
in water impact analyses and for these the properties are assumed constant at the 
following values. 

Modulus of Elasticity 

E = 29.5 x 10 psi 

Poisson's Ratio ' 

yoc = 0.30 

Specific Weight 

f = 0.283 lb/in. 3 

Also, in examination of the pertinence of core shear deformation effects: 

TP T 

Hi — _L<J jJDl 

4 

G c = 2.8 x 10 psi (Vertical shear -both directions). 

For thermal stress analyses, where the significance of the temperature dependence of 
material properties must be taken into account, the data represented in Figure n-4 
was employed. 

D. EFFECTS OF GRIDWORK REFINEMENT 

In order to determine the influence of gridwork refinement, a comparison of 
flexibility coefficients using two different grids was carried out. 

In defining the basic form of the grids which are to be used in the discrete element 
analyses of the heat shield, the shape of the sculpturing of the skins must be taken into 
account. Although all sculpturing contours are based upon radii drawn from a reference 
point lying outside the heat shield, a grid which is principally polar with the reference 
point as origin appears to be the most attractive from the standpoint of the overall 
planform. A grid using this concept is presented in Figure II- 5. In this case, there 
are 100 node points, each possessing five degrees of freedom (rotational freedom about 
a normal to the shell surface is basically deleted). When the boundary conditions are 
introduced, the total number of degrees of freedom is reduced to slightly below the 
maximum of 492 allowable in the Bell General Purpose Analysis Program. This grid 
then is typical of the maximum refinement in any particular area would necessitate an 
increase in the element size elsewhere. 


Report No. 7218-933004 


n-8 


<X (10 In./In./°F) Ablator E (10 PSI) Ablator 
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Figure n-4. Temperature Dependence of Material Properties 
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Figure n-5. Heat Shield Gridwork With Scalloped Face Thickness 
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TABLE H-2 

NODE POINT COORDINATES - REFINED GRID 


Node Pt. 

X 

Y 

Z 

1 

. -72.21 

0 

160.066 

2 

-65.155 

0 

163.065 

3 

-58.10 

0 

165.710 

4 

-53.407 

0 

167.282 

5 

-47.982 

0 

168.919 

6 

-38.624 

0 

171.300 

7 

-29.266 

0 

173.144 

8 

-19.908 

0 

174.468 

9 

-10.55 

0 

175.283 

10 

0 

0 

175.6 

11 

9.57 

0 

175.339 

12 

19.87 

0 

174.473 

13 

30.17 

0 

172.989 

14 

44.135 

0 

169.964 

15 

58.10 

0 

165.710 

16 

72.21 

0 

160.066 

17 

-72.167 

2.504 

160.066 

18 

-65.107 

2.504 

163.065 

19 

-58.046 

2.504 

165.710 

20 

-53.625 

3.126 

167.182 

21 

-48.253 

3.883 

168.796 

22 

r38.986 

5.184 

171.140 

23 

-29.718 

6.479 

174.114 

24 

-20.451 

7.781 

174.231 

25 

-11.184 

9.083 

175.008 

26 

-0.736 

10.546 

175.282 

27 

8.741 

11.877 

174.980 

28 

18.941 

13.309 

174.067 

29 

29.141 

14.741 

172.537 

30 

42.971 

16.678 

169.442 

31 

55.115 

18.383 

165.710 

32 

69.278 

20.368 

160.066 

33 

-72.032 

5.063 

160.066 

34 

-64.958 

5.063 

163.065 

35 

-57.879 

5.063 

165.710 

36 

-54.244 

6.083 

166.901 

37 

-49.021 

7.550 

168.450 

38 

-40.011 

10.078 

170.684 

39 

-31.001 

12.606 

172.382 

40 

-21.991 

15.134 

173.560 

41 

-12.981 

17.663 

174.227 
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TABLE n-2 (CONT) 


Coordinates 


Node Pt. 

X 

Y 

Z 

42 

-2.823 

20.512 

174.375 

43 

6.392 

23.094 

173.957 

44 

16.309 

25.877 

172.915 

45 

26.226 

28.659 

171.250 

46 

39.672 

32.431 

167.958 

47 

46.797 

34.433 

165.710 

48 

61.121 

38.451 

160.066 

49 

-71.796 

7.719 

160.066 

50 

-64.696 

7.719 

163.065 

51 

-57.585 

7.719 

165.710 

52 

-55.171 

8.738 

166.480 

53 

-50.171 

10.843 

167.931 

54 

-41.547 

14.476 

170.000 

55 

-32.922 

18.106 

171.533 

56 

-24.297 

21.737 

172.549 

57 

-15.672 

25.368 

173.050 

58 

-5.948 

29.460 

173.009 

5S 

2.872 

33.174 

172.415 

60 

12.365 

37.170 

171.176 

61 

21.858 

41.167 

169.301 

32 

34.734 

46.574 

165.710 

63 

49.343 

52.721 

160.066 

64 

-72.210 

11.973 

160.066 

65 

-64.045 

11.973 

163.065 

66 

-56.853 

11.973 

165.710 

67 

-52.259 

14.859 

166.984 

68 

-44.333 

19.833 

168.751 

69 

-36.407 

24.809 

169.985 

70 

-28.481 

29.784 

170.696 

71 

-20.555 

34.759 

170.894 

72 

-11.619 

40.366 

170.502 

73 

-3.514 

45.456 

169.578 

74 

5.212 

50.928 

167.972 

75 

13.935 

56.404 

165.710 

76 

29.260 

66.016 

160.066 

77 

-68.961 

21.416 

160.066 

78 

- 62.597 

18.140 

163.065 

79 

-55.050 

18.577 

165.710 

80 

-48.060 

24.800 

167.064 

81 

-41.068 

31.019 

167.889 

82 

-34.077 

37.240 

168.188 
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TABLE II-2 (CONT) 


Coordinates 


Node Pt. 

X 

Y 

Z 

83 

-27.086 

43.461 

167.967 

84 

-19.205 

50.474 

167.089 

85 

-12.055 

56.836 

165.710 

86 

0 

58.10 

165.710 

87 

0 

72.21 

160.066 

88 

-67.158 

26.535 

160.066 

89 

-60.776 

23.483 

163.065 

90 

-63.136 

35.045 

160.066 

91 

-51.029 

27.779 

165.710 

92 

-57.962 

43.066 

160.066 

93 

-45.854 

35.680 

165.710 

94 

-51.634 

50.480 

160.066 

95 

- 39.527 

42.582 

165.710 

96 

-44.153 

57.139 

160.066 

97 

-32.046 

48.463 

165.710 

98 

-34.336 

63.524 

160.066 

99 

-22.228 

53.680 

165.710 

100 

-24.163 

68.047 

1G0.0GG 
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The overall effects of grid refinement can, therefore, be studied by simplifying 
the mesh of Figure n-5 . One scheme is to omit alternate grid and radial lines. It 
is important, in the assessment of convergence characteristics, to utilize a regular 
and consistent refinement of the grid. This new grid (see Figure H-2) then represents 
effectively the coarsest acceptable mesh. Any further reduction in the number of 
node points would introduce elements of unreasonable size. 


Using these two grids, and with constant face thicknesses for simplicity, static 
analyses were performed and inter alia the influence coefficient matrix (inverse of the 
stiffness matrix) obtained. The direct flexibilities for all node points on the axis of 
symmetry have been selected as realistic characteristics of the structure and are 
presented in Figure n- 6 . The curves are not simple deflections under a fixed 
loading system, but are the deflections at each individual node point for a unit load 
applied at that point. 


As can be seen, the agreement between the two curves is very good, especially 
in the region in which impact loads will occur, i.e. to the right of center as drawn. 

The larger discrepancies between the curves, which occur to the left of center, are 
associated with the relatively larger discrete elements used in that region. In an 
actual problem, where the loading does not extend over the entire heat shield, the 
influence of the above differences is of reduced significance. From this comparison, 
it can be concluded that the coarse gridwork used in the dynamic and instability analysis 
investigations, described in later sections, yield acceptably accurate results. The 
more refined grid is retained in the present static design analyses in order to provide 
more detailed information about the peak stresses. 


E. ANALYSES IN SUPPORT OF DESIGN CHANGES 


Work on the subject contract was initiated during September, 1964. During the 
period between then and March 1, 1965, a succession of analyses were performed in 
attempts to resolve critical problems which arose as a series of new design conditions, 
proportions, and arrangement were to meet requirements imposed by water impact. 
Included were analyses which were solely intended to verify the discrete element ap- 
proach as an applicable tool for heat shield analysis through comparison with heat 
shield test data. 

Table n-3 summarizes the analyses performed prior to March 1, 1965. A total 
of sixteen analyses were listed. Each is identified by such factors as the gridwork 
employed, the nature and intensity of the imposed loading condition, the support con- 
ditions, etc. Results for these analyses were transmitted to NASA MSC, accompanied 
by sketches which describe the idealization, etc. Reference can be made to these 
sketches for further identification of the listed analyses . 
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TABL 


ANALYSES PRIOR 


Run 

Grid 

Pressure psi 


Impact Center 


1 

i 

50 

43.0 

z = 36 


2 

i 

100 

43.0 

z = 36 


3 

ii 

Variable 




4 

m 

Variable 





IV 





5 

V 

100 

13 

10° z 30 


6 

V 

100 

13 

15° z 44 


15 

V 

100 

13 

15° z 44 


7 

VI 

100 

13.2 

15° z 44 


14 

vn 

i 

100 

13.2 

15° z 44 


8 

VIII 

100 

10 

10° z 30 


9 

VIII 

100 

13 

10° z 30 


10 

vm 

100 

15 

10° z 3C 


11 

vni 

100 

10 

15° z 4* 


12 

VIII 

100 

13 

15° z 44 


13 

VIII 

100 

15 

15° z 44 


16 

IX 

100 

13 

15° z 44 

j 


♦Radius of Loaded Area 

The notation "Gpri..~ i: -’efers to tne use of an approximation 
for the flexibility of the conic section of the command module 
structure 
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(jg) BELL AEROSYSTEMS COMPANY 

s n-3 

TO MARCH 1 , 1965 


Face Sheets 

Boundary Conditions 

Group 

Case 

0.008 

Simple Support, Torus Spring + 



0.008 




0.008 

ATR 209 Torus Free 



0.008 

ATR 209 Torus Spring 



0.008 

One Quadrant, ATR 209 



0.008 

No Springs 

4 

1 

0.008 

With Springs 

4 

2 

0.008 

No Spring 

4 

B 

Scalloped 

Guide Grid, Fixed Support, Torus Free 

5 

A 

Scalloped 

Refined Grid, Fixed Support, Torus Free 

5 

B 

Scalloped 

Simple Support, Torus Spring 

6 

B-l 

Scalloped 

Simple Support, Torus Spring 

6 

B-2 

Scalloped 

Simple Support, Torus Spring 

6 

B-3 

Scalloped 

Simple Support, Torus Spring 

6 

B-4 

Scalloped 

Simple Support, Torus Spring 

6 

B-5 

Scalloped 

Simple Support, Torus Spring 

6 

B-6 

0.008 Const 

1 

1 

Simple Support, Torus Free 

7 




H-15 




© 

00 



qi/saqout - (g_Ol) £«n c H xa I.i 
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Figure n-6. Effects of Grid Refinement 


F. ANALYSES FOR CRITICAL IMPACT CASES 


Two general groups of analyses were performed of water impact conditions 
for the Apollo aft heat shield. One group pertains to drop tests performed by NASA 
of the BP28 item. These tests are known as tests No. 91, 92 and 93. The second 
group of analyses, known as the "water impact series", were conducted to determine 
the critical combinations of impact angle and wetted area in a general manner. 

Consider first the drop test analyses . The three tests correspond to water 
impact angles of 29.2°, 19.6° and 16.5° respectively. In order to provide information 
over the complete water impact, five circular wetted areas were considered for 
each of the impact angles. The radii were 10 in., 15 in., 20 in., 25 in. and 40 in. 

The fifteen design cases were then analyzed using the grid pattern of Figure n-5, 
with a nominal pressure of 100 psi in each case. The actual pressures for each case 
are taken from Figure n-7, which presents the pressure- time history for an impact 
velocity of 35 ft/sec. The corresponding radii of contact areas have been obtained 
from other data from MSC. The results from the analyses are then scaled by the 
appropriate factors. 

Results, in the form of deflection profiles along the axis of symmetry, are 
presented in Figures n-8 to 11-13 for all fifteen cases considered. The complete 
results, in the form of displacement components at all of the node points and stresses 
within the respective elements are listed in the computer analysis printout, furnished 
with this report. 

In the water impact series of analyses the objective was to determine the 
critical combinations of impact angle and wetted areas for the heat shield. The 
criteria which were to be applied are those of maximum displacements, maximum 
membrane stresses in the faces, and maximum shear stresses in the sandwich 
core. 


In order to cover the widest range of conditions, while keeping the amount of 
data as low as possible, it was decided that four impact angles combined with six 
wetted areas would provide a sufficient number of analyses. The impact angles used 
were 10°, 15°, 20° and 25°, and the wetted radii 10 in., 13 in., 15 in., 18 in., 20 in., 
and 40 in. The velocity of impact has been taken as 35 ft/sec. In general, the 
refined gridwork shown in Figure II- 5 was used. 

Figures 11-14 to 11-17 show the variation of deflection along the axis of 
symmetry with wetted radius and angle of impact. These show that the maximum 
deflection occurs in the region between 20 in. and 40 in. radius. 

The maximum stress in the faces is shown in Figure n-18. These maximum 
stresses are obtained as the sum of membrane and bending stresses and occur at 
various points on the center line of the heat shield. 
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Deflection 



Figure 11-11. Deflection Along Load Axis of Symmetry 4 >= 19.6 











Deflection 












Deflection 
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Figure n~15 # Deflection Along Load Axis of Symmetry 4 * = 15 



Deflection - In. 
























At the bolt circle, where the largest shear stresses are recorded, two shear 
stresses can be computed. These correspond to the heat shield inside and outside 
the bolt circle respectively. Since the discrete element method generates shear 
forces at the centroids of the elements, it is necessary to perform an interpolation- 
extrapolation on the element shear forces and the edge reactions to obtain the shear 
force distributions at the bolt circle itself. (The edge reaction force is equal to the 
change in shear force of the bolt circle.) 

Figures n-19, n-20, and n-21 show some typical shear stress distributions 
around the bolt circle. As indicated, the distance between the curves represents the 
bolt circle reaction. 

Figure n-22 depicts the maximum shear stresses inside and outside the bolt 
circle respectively for the range of impact centers and radii considered. 

Figure 11-23 and Figure 11-24 present the variation of the maximum face stresses 
and shear stresses with impact angle. 

For cases where the load center lies within or on the bolt circle, 20°) 
the maximum stress is 112,000 psi and lies well within the maximum allowable 
stress of 150,000 psi. For cases where the load center lies outside the bolt circle 
(|^>20 o ) the maximum stress is 147,000 psi. This can be regarded as pessimistic 
value, however, since the torus of the heat shield is assumed free to displace in 
the direction of the applied load. The effect of the actual support is to diminish 
the above predicted stress level. 

In the case of the shear forces , the peak values for both inside and outside the 
bolt circle occur when the impact center lies on the bolt circle. The maximum values 
are 3050 lb/in. and 3200 lb/in. respectively. Both values, unfortunately, lie above 
the quoted allowable of 2600 lb/in., but test results apparently indicate that the core 
can sustain such shears . 
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Figure IE-22. Maximum Core Shear Inside and Outside the Bolt Circle 
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G. THERMAL STRESS ANALYSES 


Two thermal stress analyses were performed of an idealization of the heat 
shield in which it is assumed to possess constant face sheets of 0.008 in. The heat 
shield is assumed to be hinged at the bolt circle, thereby producing a restraint of 
displacement in the radial direction of the plane of the bolt circle. The gridwork for 
these analyses appears in Figure n-25. Note that the torus section ("overhang") is 
included in the analytical representation and that the outer periphery of die torus is 
hinge supported. The temperature of the structure in the stress free state is spe- 
cified as 185°F. No external loads are considered in the two analyses described in 
this section. 

The first of the two analyses performed is designated as the "thermal soak" 
case. In this condition, a stabilized temperature of -150°F throughout the entire 
structure is realized. Since the ablator and the heat shield are composed of dis- 
similar materials, the change in temperature produces not only thermal forces but 
also thermal moments. 


The second case is designated as the thermal reentry condition. Here, the 
temperature varies over the surface of the heat shield. The outer surface of the heat 
shield. The outer surface of the ablator possesses a constant temperature of 1000° F. 
This temperature is assumed to vary across the ablator according to a parabolic law, 
reaching the same temperature as the skin of the sandwich on the other surface. The 
temperature of the sandwich is constant between the faces, but varies over the sur- 
face of the heat shield. This variation, which characterizes the temperature state, 
is defined by the listed node point temperatures, given in Table III-5. 


A determination of thermal forces and moments was necessary for both cases. 
The formulas for thermal force and thermal moment, per inch, are 


Thermal Force 


Thermal Moment 


where 


r Uo 

/To£T~ 

v 


\ 

J ^ &C / ( i? ' J 

~ o 

f JT 

T - ± 

~ c r<(s _ _ , 

1 <**■ 


The temperature (T) in these equations is the temperature change from the 
stress free state. A temperature distribution sketch follows. 
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TABLE H-4 


NODE POINT COORDINATES - THERMAL STRESS ANALYSIS GRID 


Coordinates 



X 

Y 

Z 


X 

Y 

Z 

1 

21.40 

0 

174.292 

38 

0 

71.980 

160.168 

2 

36.509 

0 

171.763 

39 

0 

75.050 

156.778 

3 

45.449 

0 

169.616 

40 

0 

74.640 

152.223 

4 

54.264 

0 

167.006 

41 

-8.189 

19.771 

174.292 

5 

58.038 

0 

165.731 

42 

-13.972 

33.730 

171.763 

6 

71.980 

0 

160.168 

43 

-17.392 

41.989 

169.616 

7 

75.050 

0 

156.778 

44 

-20.766 

50.133 

167.006 

8 

74.640 

0 

152.223 

45 

-22.210 

53.620 

165.731 

9 

19.771 

8.189 

174.292 

46 

-27.546 

66.500 

160.168 

10 

33.730 

13.971 

171.763 

47 

-28.720 

69.337 

156.778 

11 

41.989 

17.392 

169.616 

48 

-28.563 

68.958 

152.223 

12 

50.133 

20.766 

167.006 

49 

-15.132 

15.132 

174.292 

13 

53.620 

22.210 

165.731 

50 

-25.816 

25.816 

171.763 

14 

66.500 

27.546 

160.168 

51 

-32.137 

32.137 

169.616 

15 

69.337 

28.720 

156.778 

52 

-38.370 

38.370 

167.006 

16 

68.958 

28.563 

152.223 

53 

-41.039 

41.039 

165,731 

17 

15.132 

15.132 

174.292 

54 

-50.898 

50.898 

160.168 

18 

25.816 

25.816 

171.763 

55 

-53.068 

53.069 

156.778 

19 

32.137 

32.137 

169.616 

56 

-52.778 

52.778 

152.223 

20 

38.370 

38.370 

167.007 

57 

-19.771 

8.189 

174.292 

21 

41.039 

41.039 

165.731 

58 

-33.730 

13.971 

171.763 

22 

50.898 

50.898 

160.168 

59 

-41.989 

17.392 

169.616 

23 

53.068 

53.068 

156.778 

60 

-50.133 

20.766 

167.006 

24 

52.778 

52.778 

152.223 

61 

-53.620 

22.210 

165.731 

25 

8.189 

19.771 

174.292 

62 

-66.500 

27.546 

160.168 

26 

13.972 

33.730 

171.763 

63 

-69.337 

28.720 

156.778 

27 

17.392 

41.989 

169.616 

64 

-68.958 

28.563 

152.223 

28 

20.766 

50.133 

167.006 

65 

-21.400 

0 

174.292 

29 

27.210 

53.620 

165.731 

66 

-36.509 

0 

171.763 

30 

27.546 

66.500 

160.168 

67 

-45.449 

0 

169.616 

31 

28.720 

69.337 

156.778 

68 

-54.264 

0 

167.006 

32 

28.563 

68.958 

152.223 

69 

-58.036 

0 

165.731 

33 

0 

21.400 

174.292 

70 

-71.980 

0 

160.168 

34 

0 

36.509 

171.763 

71 

-75.050 

0 

156.778 

35 

0 

45.449 

169.616 

72 

-74.640 

0 

152.223 

36 

0 

54.264 

167.007 

73 

0 

0 

175.6 

37 

0 

58.036 

165.731 
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TABLE H-5 

REENTRY TEMPERATURES AND ABLATOR THICKNESSES 


Node 

Point 

Temperature 

°F 

Ablator 

Thickness in. 

Node 

Point 

Temperature 

°F 

Ablator 

Thickness 

1 

235 

1.65 

38 

162 

1.68 

2 

291 

1.54 

39 

94 

1.64 

3 

373 

1.50 

40 

80 

1.15 

4 

440 

1.49 

41 

152 

1.76 

5 

460 

1.49 

42 

142 

1.76 

6 

488 

1.43 

43 

138 

1.76 

7 

500 

1.14 

44 

128 

1.81 

8 

200 

0.80 

45 

122 

1.83 

9 

221 

1.66 

46 

110 

1.85 

10 

290 

1.56 

47 

82 

1.88 

11 

350 

1.54 

48 

80 

1.37 

12 

450 

1.53 

49 

143 

1.88 

13 

465 

1.52 

50 

129 

1.93 

14 

447 

1.46 

51 

121 

1.99 

15 

275 

1.16 

52 

112 

2.04 

16 

129 

0.84 

53 

106 

2.05 

17 

211 

1.68 

54 

96 

2.07 

18 

264 

1.58 

55 

82 

2.05 

19 

325 

1.56 

56 

80 

1.56 

20 

408 

1.55 

57 

139 

1.98 

21 

410 

1.55 

58 

120 

2.13 

22 

371 

1.50 

59 

110 

2.21 

23 

151 

1.25 

60 

101 

2.31 

24 

99 

0.90 

61 

97 

2.36 

25 

189 

1.70 

62 

90 

2.39 

26 

220 

1.62 

63 

80 

2.20 

27 

268 

1.60 

64 

80 

1.73 - 

28 

308 

1.59 

65 

137 

2.05 

29 

305 

1.60 

66 

117 

2.21 

30 

268 

1.58 

67 

107 

2i30 

31 

108 

1.45 

68 

98 

2.39 

32 

86 

1.00 

69 

93 

2.44 

33 

170 

1.72 

70 

86 

2.43 

34 

173 

1.67 

71 

80 

2.12 

35 

180 

1.65 

72 

80 

1.80 

36 

179 

1.66 

73 

185 

1.84 

37 

178 

1.67 
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ELEMENT THERMAL STRESS ANALYSIS DATA 
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Sandwich 


The dissimilarity of materials on the cross-section presents a problem in the 
determination of the proper membrane and flexural stiffnesses for the plate elements. 
Actually program input requires the plate thickness as the basic parameter with, a 
single modulus of elasticity and coefficient of thermal expansion pertinent to all 
determinations. Here, the modulus of elasticity was chosen to be constant at 
E c = 30 x 106, and the other parameters were derived as follows: 

In evaluating the foregoing integrals, the appropriate temperatures, ablator 
thicknesses, etc., are utilized for each node point. For example, the values obtained 
at the crown (node point 73) for both analyses are: 



Cold Soak 

Thermal Reentry 

Kxy 

0.94224 (10 6 ) 

0.49642 (10 6 ) 

K 

z 

0.262616 (10 9 ) 

0.108068 (10 9 ) 

F_ 

T 

-3426.24 

55.72 

M 

-1503.90 

101.06 

T 

-150. 0°F 

+185. 0°F 

T 

a 

1.84 in. 

1.84 in. 

Using these parameters, equivalent plate thicknesses are calculated. From the 
determined thermal force, an equivalent coefficient of thermal expansion (o£ c ) was 
then derived; retaining the actual temperature difference from the stress free state 
(AT). For node point 73 the equivalent values are: 


Cold Soak 

Thermal Reentry 

t 

m 

s 

oc 

c 

0.03141 

0.8179 

0.2714 (10~ 8 ) 

0.01655 

0.6084 

0.1121 (10 -8 ) 

M«c 

-1503.9 

101.06 

At 

-335°F 

1°F 
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With these equivalent section and material properties, the general purpose 
program then computes the correct stiffness and thermal forces. 

Results of the thermal analyses appear in Figure 11-26 and Table II- 7 . The 
figures illustrate the variation of displacement along the axis of symmetry. The 
specific values of node point displacements upon which these representations are 
based are given in Table II- 7 , as are the displacements at other selected points 
on the heat shield. 

Figures n-27 to 11-31 present the variation of meridional and circumferential 
thermal stresses in the heat shield along the axis of symmetry. 

In Figures n-27 to n-30 the stresses in the upper and lower faces of the 
sandwich are shown for the thermal reentry and cold soak conditions. The peak 
stress of 122000 psi occurs in the upper face near the bolt circle during the cold 
soak. This stress, although high, is acceptable. 

For the ablator the stress in the upper edge for the cold soak condition 
(Figure 11-31) reveals that the maximum stress is 1050 psi (tensile). This is well 
above the maximum strength of the ablator material (approximately 600 psi) . 

However, in evaluating these results, especially the above ablator stress, the 
entire thermal stress analysis must be carefully examined. In particular the boundary 
conditions must be considered. At the bolt circle, it has been assumed that all 
translational movement was completely restrained. In actual practice, the heat 
shield is attached to the (non-rigid) internal cabin structure by supports which also 
exhibit some degree of flexibility. While this assumption may be satisfactory for 
conditions in which the loading is principally normal to the surface of the shell, in 
the thermal cases, the radial expansive forces will probably cause appreciable 
movement at the bolt circle with consequent reduction of the stresses quoted above. 

In addition, no account has been taken in the analysis of the built-in thermal 
stresses which exist in the heat shield. These arise from the fact that although the 
heat shield zero stress temperature is taken as 185°F (the temperature at which the 
ablator is attached), the assembly of the complete heatshield to the support structure 
actually occurs at around 80° F 

The thermal stresses obtained above may be regarded as having been derived 
by a realistic approach to an unrealistic problem. In order to obtain more realistic 
results it would be necessary to include the inner cabin structure in the analysis. To 
do this would, in turn, increase the magnitude and complexity of the problem well 
beyond the limits of the current contract and cannot therefore be considered at 
present. 
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Deflection - Inches Deflection - Inches 



b. Thermal Soak 

Figure n-26. Thermal Displacements 
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TABLE DISPLACEMENT NORMAL TO SURFACE - THERMAL CASES 
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m. DYNAMIC RESPONSE 


A. INTRODUCTION 

The basic design analyses described in the previous chapter have been per- 
formed on a purely static basis, utilizing very refined gridworks. The deflection 
curves, which have been obtained for various radii of wetted areas corresponding to 
points of the water impact pressure-time history do not include dynamic effects. 
Since the times involved are very small (<0.05 sec), inertial and dynamic response 
effects may have an appreciable significance. Analyses in which the effects of 
dynamic response have been neglected generate results which are expected to be 
conservative. It is desirable to obtain some quantitative estimate of the degree of 
conservatism, if any, involved in the previous static analyses. For this purpose, 
the dynamic response of the Apollo heat shield due to water impact was determined, 
and the results are compared with those obtained from the static design analyses 
in this section. 

The items of work covered by the complete dynamic response analyses pro- 
gram were as follows: 

(1) After review of the results of the static analyses, loading cases and an 
idealization for the heat shield were defined. 

(2) Using the gridwork selected, an analysis was performed to determine 
the displacement influence coefficient matrix. A mass matrix for the 
heat shield was also derived. 

(3) For the selected analysis conditions and structural idealization, load- 
time histories, F(t), were developed for all loading points. 

(4) Dynamic analyses were then performed using the forcing functions 
defined previously. 

(5) Estimates were made of the significance of dynamic, as compared with 
static, behavior with regard to stresses. 

Each of the above items of work is discussed in the following sections. A detailed 
explanation of the theory and procedures underlying these dynamic analyses appears 
in Appendix B. 

B. STRUCTURAL IDEALIZATION 

The chosen idealization appears in Figure III-l. In this idealization, a 
total of 97 grid points on the surface (or middle surface) of the heat shield are free 
to displace. Although five degrees of freedom at such points appear in the deter- 
mination of stiffness, only one degree of freedom at each point - that which is 
normal to the surface — is considered in dynamic response analysis. 
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The procedure for obtaining the desired flexibility characteristics is as follows. 
The stiffness matrix for the structure is constructed by the Bell General Purpose 
Structural Analysis Program in the form 



where % , ^ and ^ are coordinate directions in and normal to the tangent planes at the 
respective points and 0 and 0 are angular displacements in the * and ^ directions. 
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(ni- 2 ) 


n UCl C UiC V A ai V I1UAIUHU.J 

however, are only those which relate the Px’$ to the i y in 


In the Bell General Purpose Structural Analysis Program, instructions are 
provided which allow the assignment of the complete inverted stiffness matrix to 
tape storage. The analyst then designates the rows and columns to be recalled and 
thereby constructs the desired flexibility matrix, in this case . 


The idealization shown in Figure in-1 was also used to determine the mass 
properties of the heat shield. Initially the surface area subtended by sets of grid 
points was calculated, such as the area bounded by points 11, 12, 19, and 18. One- 
quarter of this area was assigned to each corner grid point. In this manner the 
surface area assigned to each grid point is found. Next, the ablator, face sheet and 
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core thicknesses were determined for the assigned area resulting in assigned volumes 
for each of these structural elements. These volumes were multiplied by the proper 
density values to obtain weights. Table III-l lists the assigned areas, structural 
element weights and total weight. 

The mass matrix[ M3, to be used in the dynamic response study is easily con- 
structed from Table in-1 by placing the pertinent weight, divided by f , for each grid 
point on the main diagonal of this matrix. This approach, commonly labeled the 
"mass lumping technique,” results in zero-valued off diagonal terms. 

C. LOAD-TIME HISTORIES 

Forcing functions or load-time histories, F(t), were developed for each grid 
point. Each loading curve is based on the impact pressure-time and wetted radius- 
time histories shown in Figure DI-2 for a chosen impact angle of = 15° and vertical 
impact velocity of 35 ft/sec. A particular grid point load-time history was con- 
structed by determining the time , , at which the grid point was "wetted" through 

use of the wetted radius-time history curve and grid point coordinates. The pressure- 
time history curve is entered at-*f» to find the pressure decay for the grid point. Load 
time histories are easily determined by multiplication of the grid point pressure- 
time history by the pertinent surface area. Typical results are shown in Figures 
m-3 and ID-4 for grid points 7 and 57 respectively. These results are conservative 
in nature since the heat shield is to be subjected (theoretically) to a peak force which 
is greater than that obtained if the actual rate of wetted grid point area were con- 
sidered. In the latter instance a monotonic increase of force with time would occur 
starting at some and peaking at some >tf. The peak value would be less 

than that shown due to the pressure decay. This phenomenon is of negligible impor- 
tance for grid points near the impact point where the radius of wetted area changes 
rapidly with time. 

The forcing function data are prepared as input to the dynamic response pro- 
gram by providing the coefficients of the exponential function 

F(t) - ,4 e B * +C (m-3) 


These coefficients were determined by a "least-squares" technique applied to the 
basic pressure decay curve of Figure m-2. The most accurate results were obtained 
by fitting exponentials to the time regions 0.00l£ t 0.010 and 0.010 < t 0.05. 
These are shown on Figure HI-2. The coefficient A in Equation HI-3 represents the 
assigned grid point area. It is noted that Equation m-3 is only used in the time 
region . For Jtp F(t) is set equal to zero. 
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TABLE m-1 


NODE POINT COORDINATES - REVISED GRID 


(See Figure III-l for Layout) 



X 

Y 

Z 


X 

Y 

1 

-72.21 

0 

160.066 

41 

-29.012 

16.75 

2 

-65.155 

0 

163.065 

42 

-20.785 

12.0 

3 

-58.100 

0 

165.710 

43 

-11.184 

9.083 

4 

-51.00 

.0 

168.031 

44 

-58.419 

42.444 

5 

-43.0 

0 

170.254 

45 

-47.004 

34.150 

6 

-33.500 

0 

172.375 

46 

-41.260 

29.977 

7 

-24.000 

0 

173.952 

47 

-34.788 

25.275 

8 

-10.550 

0 

175.238 

48 

-27.102 

19.691 

9 

-71.814 

7.548 

160.066 

49 

-51.634 

50.480 

10 

-64.798 

6.810 

163.065 

50 

-39.527 

42.582 

11 

-57.781 

6.073 

165.710 

51 

-34.077 

37.240 

12 

-50.721 

5.331 

168.031 

52 

-28.481 

29.784 

13 

-42.764 

4.495 

170.254 

53 

-44.153 

57.139 

14 

-33.316 

3.502 

172.375 

54 

-32.046 

48.463 

15 

-23.868 

2.509 

173.952 

55 . 

-27.086 

43.461 

16 

-70.632 

15.013 

160.066 

56 

-20.555 

34.759 

17 

-63.731 

13.546 

163.065 

57 . 

-15.672 

25.368 

18 

-56.830 

12.080 

165.710 

58 

-12.981 

17.663 

19 

-49.886 

10.603 

168.031 

59 

-34.336 

65.524 

20 

. -42.060 

8.940 

170.254 

60 

-22.228 

53.680 

21 

' -32.768 

6.965 

172.375 

61 

-19.205 

50.474 

22 

-68.676 

22.314 

160.066 

62 

-11.619 

40.366 

23 

-61.966 

20.134 

163.065 

63 

-5.948 

29.460 

24 

-55.256 

17.954 

165.710 

64 

-2.823 

20.512 

25 

-48.504 

15.760 

168.031 

65 

-0.736 

10.546 

26 

-40.895 

13.288 

170.254 

66 

-24.163 

68.047 

27 

-31.860 

10.352 

172.375 

67 

-12.055 

56.836 

28 

-22.825 

7.416 

173.952 

68 

-3.514 

45.456 

29 

-10.643 

3.458 

175.243 

69 

2.872 

33.174 

30 

-65.967 

29.370 

160.066 

70 

6.392 

23.094 

31 

-59.522 

26.501 

163.065 

71 

8.741 

11.877 

32 

-53.077 

23.631 

165.710 

72 

0 

72.210 

33 

-46.591 

20.744 

168.031 

73 

0 

58.100 

34 

-39.282 

17.490 

170.254 

74 

5.212 

50.928 

35 

-30.604 

13.625 

172.375 

75 

12.365 

37.170 

36 

-62.536 

36.105 

160.066 

76 

16.309 

25.877 

37 

-56.426 

32.578 

163.065 

77 

18.941 

13.309 

38 

-50.316 

29.050 

165.710 

78 

29.260 

66.016 

39 

-44.167 

25.500 

168.031 

79 

13.935 

56.404 

40 

-37.239 

21.500 

170.254 

80 

21.858 

41.167 
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172.375 

173.952 

175.008 

160.066 

165.710 
168.031 
170.254 

172.375 

160.066 

165.710 
168.188 
170.696 

160.066 

165.710 
167.967 
170.894 
173.050 
174.227 

160.066 

165.710 
167.089 
170.502 

173.009 

174.375 
175.282 

160.066 

165.710 
169.578 
172.415 
173.957 
174.980 

160.066 

165.710 
167.972 
171.176 
172.915 

174.067 
160.066 

165.710 
169.301 


HI-5 


TABLE m-1 (CONT) 



Y Z 



81 

26.266 

28.659 

171.250 

82 

29.141 

14.471 

172.537 

83 

49.343 

52.721 

160.066 

84 

34.734 

46.574 

165.710 

85 

39.672 

32.431 

167.958 

86 

42.971 

16.678 

169.442 

87 

61.121 

38.451 

160.066 

88 

46.797 

34.433 

165.710 

89 

55.115 

18.383 

165.710 

90 

69.278 

20.368 

160.066 

91 

72.210 

0 

160.066 

92 

58.100 

0 

165.710 

93 

44.135 

0 

169.964 

94 

30.170 

0 

172.989 

95 

19.870 

0 

174.437 

96 

9.570 

0 

175.339 

97 

0 

0 

175.6 
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D. DYNAMIC RESPONSE 


The dynamic response of the Apollo heat shield was determined by use of the 
techniques described in Appendix B. This technique utilizes the natural modes and 
frequencies of vibration of the subject structure. 

A matrix iteration technique was used to determine the vibrational character- 
istics. Mathematically, the eigenvalues (natural frequencies) and eigenvectors (mode 
shapes) of the following matrix equation must be found. 

CS.J < m - 4 > 

The matrix! £**Jis the set of flexibility influence coefficients described by Equation 
(III-4)j$|^is a column of displacements and£ M 3 is the diagonal mass matrix con- 
structed from the total weight column shown in Table IE-1. The flexibility matrix is 
reduced to 80th order by virture of the boundary conditions along the centerline of 
the heat shield and bolt circle support points. The mass matrix is also of this 
order. The first ten modes and frequencies were obtained from Equation (IH-4) and 
these are tabulated in Table ID-2. Additional mode shapes are given in Reference 
. The mode shapes for the first two modes are shown in Figures III-5 and III-6 
where negative modal values at each grid point are shown by the cross-hatched 
areas. Relative amplitudes for the centerline grid points, using grid point 97 as the 
normalizing station, are shown in the upper portion of each figure. Approximate 
node line locations (points of zero displacement) are also shown. Examination of 
these figures shows that only modes symmetrical about the heat shield centerline are 
obtained. This is due to the fact that the flexibility matrix [ $ xx ] was constructed on 
the basis of symmetric loadings being applied to the subject structure. The 
associated modal deflection pattern is reflected in the boundary conditions consistent 
with the applied loadings. The calculated frequencies associated with each mode are 
shown in Table ni-2. 

TABLE III- 2 

VIBRATION MODE AND FREQUENCY DATA 



Frequency 

Generalized 

Generalized 

Generalized 

Mode 

(rad/sec) 

Mass 

Stiffness 

Damping 

1 

592.0 

19.80 

6.94 x 10* 

100.41 

2 

657.6 

47.50 

20.56 x 10 6 

297.55 

3 

683.6 

12.10 

5.67 x 10® 
28.92 x 10 6 

82.04 

4 

750.6 

51.30 

418.45 

5 

813.9 

187.00 

123.84 x 10 6 

1719.94 

6 

848.7 

0.10 

0.07 x 10 6 

1073.55 

7 

869.9 

4.30 

3.28 x 10 6 

47.49 

8 

941.6 

42.17 

37.39 x 10 6 

540.96 

9 

984.1 

0.68 

0.66 x 10 6 

9.54 

10 

1112.0 

128.39 

158.76 x 10 6 

227.29 
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Figure HI-5. Heat Shield First Vibration Mode Shape (fj^ = 94.26 cps) 



Figure HI-6. Heat Shield Second Vibration Mode Shape (f 2 = 104.74 cps) 


The frequencies range from 592 rad/sec (94.3 cps) to 1112 rad/sec (178.5 cps). 


The method of solution for dynamic response as described in Appendix B 
assumes that the displacement of any grid point due to a dynamic load can be 
obtained by a summation of products of mode shapes (relative displacements) and 
generalized coordinates. That is 


j ^ ^ i ^ ■> <f‘ ^ 


(m-5) 


where y£C ( x c - ft ) is the jth normalized vibration mode shape evaluated at the 
ith point and f^.<«)is the generalized coordinate in the jth mode. Use of Equation 
(D3-5) leads directly to generalized mass, stiffness and damping terms defined as 
follows 



I 




/*? , 


l ( X,; 



(m-6) 


ICj,= cjj Mj <ni-7) 



In the present application the summations extend from i = 1 to 80, j = 1 to 10. 
Generalized mass, stiffness, and damping values associated with each mode are 
listed in Table m-2. These values form the elements of the diagonal mass, stiffness 
and damping matrices in the equations of motion. (See Equation B-23 , Appendix B.) 

Although formation of the generalized mass and stiffness terms is straightfor- 
ward, some discussion of the damping terms is needed. The coefficient f'i in 
Equation (in-8) is defined to be the damping coefficient associated with each vibration 
mode and is expressed in terms of critical damping as 


fi ■■ where C c 1 2 Vk Mj (HI-9) 

In the present application was assumed to be equal for all modes. The parameter 
CO in Equation (HI-8) is the frequency response of each generalized coordinate and is 
unknown at the outset of the analysis. This frequency is easily estimated however 
by averaging the response frequencies for the undamped case. For the present 
analyses, CO was assumed to be 110.0 cps. 
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The matrix equations of motion (Equations B-23 , Appendix B) were solved by 
the use of the Bell Matrix Iteration Time History Program. Program output consists 
of time histories of the generalized coordinates and total grid point displacements. 
Typical results are shown in Figures III-7 to III-10 (additional results are displayed 
in Reference 2 ). The generalized responses shown are to be expected since the 
ratio of the period of the generalized force to natural period ranges from five to nine 
approximately. Upon removal of the exciting forces at t = 0.050 seconds each 
generalized mass-spring system oscillates at its natural frequency. The effect of 
structural damping upon maximum displacement is negligible for g = 0.01 as shown 
by the open circular symbols on Figures ni-7 and ni-8. An extreme case of damp- 
ing, g = 0.40, shows a substantial decrease in amplitude and for £ g the motion 
is entirely damped-out before the excitation pulse has vanished. 

Total grid point displacements for zero damping are shown in Figures (III— 9) 
for the centerline grid points 4 and 97. These are typical of the results obtained and 
additional grid point time histories are presented in Reference 5. Examination of 
the results showed that the frequency response during excitation varied between 100 
and 140 cps approximately. During the period 0.05 — 0.080, when the external 

forces are assumed to be zero, the frequency response in approximately 125 cps for 
all grid points. 

A comparison between static and dynamic response is shown in Figures ni-10 
and III— 11 at selected time intervals for the centerline grid points. These figures 
show that the deflections obtained by static response, in most instances, is greater 
than that obtained by dynamic response analyses. Of particular interest is the fact 
that as time increases the dynamic and static responses on the lower half of the heat 
shield (negative coordinate) become closer in agreement. Thus, the dynamic effects 
represent stress levels lower than the levels predicted in a static analysis. This 
conclusion lends justification to the performance of static analyses and loadings of 
the maximum acceleration, rather than for the apparently more critical loadings 
earlier on impact. 
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Figure III-8. 
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Figure HI-10. Displacement Along Load Axis of Symmetry 
Impact Angle <t> = 15° 















- Dynamic Response (g = 0) 
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IV. INSTABILITY ANALYSES 


A. OBJECTIVES 

The water impact condition produces membrane stresses throughout the heat 
shield and it is conceivable that these stresses could be large enough to produce a 
failure of the heat shield in an elastic instability mode. The objectives of the work 
described in this section were to obtain the critical values of applied load for the 
water impact design condition based on a discrete element approach to linear 
instability analysis, modified by empirical correction factors. 

The problem of determining elastic instability effects on the behavior of the 
scalloped heat shield during the water impact phase is one which cannot be solved 
through use of classical techniques with the expenditure of modest amounts of time 
and effort. Developments of the matrix displacement method, however, permit an 
engineering solution of this problem within the framework of this study. 

The detailed work items in this area were as follows: 

1. Perform instability analyses of a non-scalloped heat shield configuration 
for fixed support at the bolt circle under a uniform load condition. 

2. Define available classical solution to the above case. 

3. Compare results of discrete element and classical solutions. 

4. Perform analyses of scalloped heat shield for water impact design 
conditions. 

Item (4) was conducted for the case of a 10° impact angle and wetted radii of 10, 
20, and 40 inches, respectively. A description of these efforts is given in the follow- 
ing paragraphs. 

B. DESCRIPTION OF PROCEDURE 

A detailed development of the present approach to elastic instability analysis 
is given in Appendix C . The following is a review of the basic concepts of this 
approach and its illustration by means of a simple problem. 

It is possible to define, in connection with elastic instability phenomena, a 
hierarchy of sophistication in achieving solutions. This hierarchy is as follows: 
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a. Frame 


b. Idealization 


Figure IV- 1. Simple Frame for Illustration of 
Instability Analysis Procedure 
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(1) The middle surface direct stresses are evaluated in a membrane analysis. 
This approach is common to "classical" methods of instability analysis. 

It is not adopted in the discrete element approach because: 

(a) The middle surface direct forces are as conveniently determined in 
an analysis which includes flexure as in one which does not. 

(b) It is possible for a singular stiffness matrix to be associated with a 
purely membrane model of a structure which is designed to carry 
bending moments . 

(c) It is not convenient to input predetermined values of the internal 
forces into the existing Bell buckling analysis program. 

(d) In certain cases, the membrane stresses are significantly affected 
by flexure. 

(2) The middle surface stresses are determined in an analysis that includes 
both flexure and axial behavior, but excludes finite deflection effects. 

(The axial strain-displacement equations do not include strains due to 
flexure.) 

(3) The determination of middle surface stresses stems from an analysis 
which includes finite deflection effects . 


Procedure (2) is perhaps clarified by the problem shown in Figure IV -la, 
which is a simple frame. It can also be regarded, however, as the two-dimensional 
counterpart of the shell problem. Figure IV-lb shows the discrete element idealiza- 
tion. Only point 2 is free to move so that stiffness equations need be written only for 
that point. Since instability effects are at first disregarded, in determining the axial 
force distribution, we can write these equations as follows: 


F* * 

r X 

Q 

L L 

F»* ' 

■K -- 

M w, 

f* = *>• - e 


-C.EI + e 

n“= - 4 -S. w , + - !tO 

¥ 

L L 

7 L*- L 



( xsz: - 1) 


Report No 


7218-933004 


IV- 3 


This may not always be the case, however. Furthermore, the above simple 
problem represents the hemisphere, rather than the shallow spherical shell. To 
examine this more closely, a series of analyses was conducted for the structure 
shown in Figure IV-2 (with the same cross-sectional properties as before) and the 
accompanying table was constructed to indicate the % improvement over axial loads 
computed from "membrane theory" 


P = 1000 lb 



Figure IV-2. Frame for Numerical Analysis 


TABLE IV- 1 

EXAMPLE FRAME RESULTS 


Degrees 

Uncoupled 

Membrane 

Flexure 

Behavior 

Coupled 

Membrane 

Flexure 

Behavior 

% 

Difference 

45 

707.0 

702.2 

0.8% 

30 

iooo:o 

980.0 

2.0% 

15 

1932.0 

1761.0 

8.9% 

5 

5737.0 

3008.0 

47.5% 
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where 


cnr-o 


Let us examine the improvement in the predicted element membrane forces 
which is realized by utilization of the above formulation rather than the "membrane" 
theory which is represented by 


P- ^ 
F - + 


(XSC-fo) 


Consider the case defined by 

? K = 1000 lb 

P* = 0 

My= 0 

4 = 45° 

0 

From Equation IV -6 we have Fj, - 1000 lb. Alternatively, using the previously 
derived equations, we have 

6 A£ 

c - — ui 

' X - L 

= C -fe) 


(m -7) 


i A -») 


" 1,er * X = ( t + 

With L - 12 inches and with a square cross section, lxl, the result is 


= <ns\i 


Hence, the difference between this and the membrane solution is not significant. (0.5%) 
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Interpretation of these results for the elementary frame structure enables 
general conclusions to be drawn regarding the nature of discrete element idealiza- 
tions of practical structures. Of primary interest is the demonstrated need for the 
consideration of membrane-flexure coupling. The importance of this consideration 
for shallow configurations of the example structure is apparent from the percentage 
error data shown in Table IV-1. Accordingly, interelement coupling of membrane 
and flexure behavior was incorporated in the instability analyses conducted for the 
Apollo heat shield structure. A detailed description of the instability analysis method 
developed is presented in Appendix C. 

C. EVALUATION ANALYSES 

Two classical problems have been solved with the subject discrete element 
approach to elastic instability analysis in order to verify its accuracy. One of the 
problems involves the prediction of the buckling pressure for a uniformly loaded 
circular arch. The other is concerned with the prediction of the deflections and the 
linear buckling of a pressurized spherical cap. 

The first problem is illustrated in Figure IV -3. A circular arch, subtending 
an angle at 120°, is subjected to a uniform normal pressure (p, lb/in.). The arch is 
fully fixed at the ends and possesses a cross-sectional area A, a moment of inertia I, 
and is made at a material with a modulus of elasticity, E. 

The arch was idealized as a system of 12 beam-column elements (Figure IV-3). 
The eigenvalue was computed to be 

o * 8 El A 

pcr IT 

while the corresponding buckling mode was of the nonsymmetric form shown in 
Figure IV -3b. 

This problem was solved by Wempner and Kesti, using classical techniques, 
in Reference 3 . They obtain the solution 

pc - *** 
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a. Problem Data 



Figure IV-3. Deep Ar 


which differs from the computed solution by less than 1%. Also, Wempner and Kesti 
predict a mode shape which is essentially identical with that shown in Figure IV-3b. 

The close agreement between the numerical and theoretical results in this case 
indicates that the linearized discrete element method used here is indeed applicable 
to this class of problem. In this case, where the subtended angle is relatively large, 
the effects of large deflections (as characterized by terms of the type 1/2 
have only a secondary effect, and thus the linear theory is applicable. In very shallow 
arches, these large deflection terms, can assume dominant importance and an anal- 
ysis disregarding such effects may yield a very poor approximation. 

In order to check the validity and applicability of the discrete element method 
in the prediction of the buckling behavior of shallow spherical caps similar to the 
Apollo heat shield, experimental results obtained by Kaplan and Fung (Reference 4) 
were examined. They conducted a series of tests on shallow spherical caps subjected 
to uniform hydraulic or pneumatic pressure. In these tests, not only the critical 
pressures were obtained, but also the deflected shapes were measured in a number 
of cases up to collapse. 

Their Specimen 21 was employed for the purposes of comparison. The 
principal dimensions and properties are indicated in Figure IV -4b. The cap is fully 
fixed at the edge and subject to air pressure. 

In the numerical analysis, only one quadrant of the cap is considered. The 
gridwork used is shown in Figure IV-4b. For this shell, the critical pressure 
has been determined experimentally to be 34.40 psi. The deflected shape has been 
measured and is presented in Reference 4 for a pressure of 0.581 ft r . Thus, in 
the numerical analysis the loading has been taken as 

0.581 * •* 

Using this pressure, the deflected form shown in Figure IV-4b was determined. 
The comparison between the experimental and theoretical results is excellent. In 
the same analysis an apparent buckling pressure of 63.1 psi was determined by an 
eigenvalue routine. Classical linear shell buckling theory yields a critical pressure 
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0.008 in. 



Figure IV-5. Heat Shield With Constant Face Thickness Uniform Pressure 






Gerard and Becker (Reference 5) have examined the comparison between 
experimental results and classical linear theory for this class of problem. Based 
upon the available test data, it is demonstrated that, due to the effects of initial 
imperfections and large deflections, the classical factor of 1.2 must be reduced. 

The amount of the reduction is dependent upon the R/t ratio, and for the present 
case of R/t = 362, the reduction factor is approximately 0.42. When the case of the 
sandwich heat shield is considered subsequently, a similar reduction factor is applied 
to the values obtained by the linearized numerical analysis to account for the effects 
of large deflections and initial imperfections . 

D. ANALYSIS RESULTS 

1. Idealized Heat Shield - Uniform Pressure 

The first instability analysis, which was performed upon a structure 
representative of the actual Apollo heat shield, was intended to determine the 
critical uniformly distributed pressure which causes buckling of a spherical sandwich 
cap. The cap has the same dimensions, etc., as the heat shield except that the sand- 
wich face sheets have a uniform thickness of 0.08 in. For the purposes of this 
analysis only one quadrant of the cap was considered. Full details of the cap and the 
gridwork are given in Figure IV-5. 

The use of the matrix instability analysis concepts described in Appendix 
C take into account the simultaneous action of membrane and flexural stress sys- 
tems, but the effects of flexural deformations on the element membrane strain- 
displacement equations are neglected. If then, in a problem, the membrane stress 
state undergoes a simple linear change in magnitude, but is not redistributed as the 
external applied load is varied linearly, then the critical buckling load of a system 
may be estimated reliably (ignoring large deflection terms) by factoring any applied 
load by an eigenvalue determined for that load level. 

In the case of the spherical cap, the bending deformations cause redistri- 
butions of the membrane stresses. The redistributed membrane stresses alter the 
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bending stiffness, and thus it is necessary to iterate at each load level to determine 
the actual stress distributions. In order to find the actual value of the buckling load, 
it is therefore necessary to reach the situation in which the eigenvalue is unity; 
i.e., no extrapolation is required. The applied load is then the buckling load. 

This condition may be reached by simple trial-and-error, but a more 
sophisticated approach is desirable. For the present case, a simple graphical tech- 
nique was used initially, and based upon these results a mathematical approach was 
developed. 

When the spherical cap was subjected to a uniform pressure of 10 psi, the 
apparent buckling load was computed as 84 psi (see Figure IV-5). Application of 
20 psi yielded 81 psi, etc. The curve in Figure IV-5 shows the various applied loads. 
The intersection of the curve and the line X = 1 indicates the critical value of the 
applied pressure as 47.0 psi. As discussed previously, this value must be reduced 
by a suitable factor to take account of large deflection effects. Curves giving this 
factor, presented in Reference 9 are based upon the radius-to-thickness ratio. For 
a sandwich panel, however, the definition of the thickness presents a problem, since 
face, core, or effective bending thicknesses may be used. Fortunately, in the region 
covered by these parameters, the curve of the reduction factor of Reference 4 is 
relatively flat and a mean value of approximately 50% may be selected. The critical 
uniform pressure for the idealized heat shield then becomes 

f^cr' = 0.5 .0 23.5 p.vv, 

2. Actual Heat Shield - Water Impact Cases 

In view of the iterative nature of the instability analysis it was decided 
to consider initially only three actual water impact cases. These are all for an 
impact angle of <|> = 10° and for the three radii of wetted area: r = 10”, 20", and 40". 
As discussed previously, the instability analyses utilize a considerably coarser 
gridwork than that used for the basic design analyses. The modified gridwork, which 
is shown in Figure IV-6, is based upon the polar grid of the design analyses 
(c.f., Figure III-5) and has approximately 150 degrees of freedom - the limit of the 
current instability program. v 
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Figures IV-7 to IV-9 present the results of the three instability analyses. 
The buckling pressures determined from these curves are as follows: 

Radius 10" 20" 40" 

Pressure (psi) 1179.0 396 210 

In the previous section, the concept of a reduction factor which accounts 
for the large deflection effects not included in the above analyses was introduced. 
The precise value of the factor is questionable but from data given in Reference 
it is estimated to be 0.5. From MSC data for the case of 35 ft/sec water impact 
velocity the pressures experienced at the appropriate radii can be compared with 
the above values after factoring by 0.5. 



10" 

20" 

40" 

Factored Pressures 

589.5 

198.0 

105.0 

Estimated Pressure 

273.0 

151.0 

80.0 

(MSC) 




Safety Factor 

2.16 

1.31 

1.31 


From these figures it may be seen that an overall instability failure of the Apollo 
heatshield is not predicted. 

E. FINITE DISPLACEMENT CONSIDERATIONS 

Extension of the instability analysis to include finite deformation and finite 
displacement effects could be expected to improve significantly the prediction of 
shell behavior for certain shallow configuration and applied loading situations. Such 
an extension could be effected within the framework of the instability analysis pro- 
cedure developed herein. 

The behavior of a structure idealized as an assembly of discrete elements is 
characterized by the nature of the discrete element representations. Accordingly, 
in the development of a nonlinear formulation, attention may be focused on the 
derivation of appropriately nonlinear discrete element representations. 

A nonlinear mathematical representation for a structural discrete element is 
most easily derived using the principle of stationary potential energy in conjunction 
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Critical Pressure - psi 



0 1000 2000 3000 

Applied Pressure - psi 


Figure IV-7. Heat Shield with Scalloped Face Thickness, 10 in. Radius Loaded 

Area, f = 10° 
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Critical Pressure - psi 


Loaded Area 



Applied Pressure - psi 

Figure IV-8. Heat Shield with Scalloped Face Thickness, 20 in. Radius Loaded 

Area $ = 10° 
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Critical Pressure - psi 


<£= 10°, R = 40 


Loaded Area 




Figure IV-9. Heat Shield With Scalloped Face Thickness, 40 in. Radius 

Loaded Area, (p = 10° 
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with the Rayleigh-Ritz procedure. In order to account for finite deformation and 
finite displacement effects, it is necessary to retain selected nonlinear terms in the 
element strain-displacement relations. For thin shell and slender prismatic struc- 
tural components it suffices to include second order terms in the transverse dis- 
placement functions. The retention of quadratic terms in the strain-displacement 
relations results in a quartic potential energy functional. Construction of displace- 
ment functions ih accordance with the Rayleigh-Ritz procedure and execution of the 
variation of the quartic potential energy yields a discrete element representation in 
the form of a set of nonlinear algebraic equations which may be written as: 


where 


t-M 

LU 

'■'-i 

M 


I nonlinear algebraic equations wmcn may be written as : 

+ [XHM + 

is the vector of forces applied to the element at the node points, 
is the corresponding vector of node point displacements, 
is the familiar element stiffness matrix. 

is defined as the "first order element incremental stiffness matrix" . 

Its elements are linear functions of the element node point displace- 
ments. 

is defined as the "second order element incremental stiffness matrix". 
Its elements are quadratic functions of the element node point displace- 
ments. 


It is instructive to exhibit the incremental stiffness matrices for a frame 
discrete element appropriate to the elementary frame structure considered 
previously (Figure IV-2). 
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The approach to incorporating finite deformation and finite displacement effects 
outlined above is readily applicable to the thin shell elements appropriate to the heat 
shield analysis as well as the frame element considered here. A further important 
feature of this method is that the numerical algorithm required for the prediction of 
behavior is conceptually unchanges though revisions might be necessary to obtain 
favorable convergence characteristics. 

In a review of the work described earlier, Budiansky and Sanders (Reference 18), 
strongly urge the performance of nonlinear deflection analyses and present a theoret- 
ical extension of discrete element concepts to accomplish this end. Their suggestions 
parallel the above procedures to a significant extent. 
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V. DISCRETE ELEMENT PROPERTIES 


A. HISTORICAL BACKGROUND 

In order to place this portion of the study in proper perspective, a brief 
historical review is given of the development of the discrete element approach to 
structural analysis. 

The concept of physically idealizing complex structures as assemblies of 
discrete structural elements has been known and utilized for over 100 years, in the 
form of procedures for truss and frame analysis. The potential value of this approach 
as a means for the analysis of structural continua such as plates and shells remained 
unrecognized until the early 1950's. Then, a large number of papers, were published 
which clarified and refined the concept of discrete element structural idealization 
and generalized the analysis procedures associated with such concepts. A landmark 
paper in this connection is that by Turner, et al (Reference 6). Numerous develop- 
ment efforts were initiated to construct discrete element formulations suitable for 
realistic idealization of a wide variety of practical structures. 

Many mathematical models were derived for a number of discrete structural 
elements using several distinctly different approaches. The single aspect common 
to the various investigations was a general inattention to the convergence character- 
istics of the discrete element representations put forward. Most element repre- 
sentations were found to exhibit satisfactory convergence while others proved use- 
less in the absence of favorable convergence characteristics. Through a heuristic 
process, general analysis capabilities of broad applicability were developed by 
building up computer program libraries of the best discrete elements. 

In 1963, Melosh (Reference 7), laid down guidelines, based on established 
principles of structural mechanics, for achieving convergent discrete element 
representations. The definition of procedures for developing element representations 
within the confines of these guidelines proved to be an elusive goal. As a consequence, 
the huge strides made in the development of additional and improved discrete element 
representations continued to stem from the heuristic process. 
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Only recently, in October of 1965, several investigators (8, 9, 10) have presented 
pre limin ary reports defining general systematic procedures for deriving convergent 
representations for a broad class of discrete structural elements. Those which are 
pertinent to heat shield analysis will be discussed subsequently. 

B. SANDWICH ELEMENT 

Prior to the initiation of this study, stiffness equations for the triangular 
sandwich element (Figure V-l) had been formulated at Bell Aerosystems Company. 
This element consists simply of upper and lower faces which are the conventional 
triangle in plane stress (Reference 6), separated by a core which carries only vertical 
shear stress. By combining these components to form a single element, including 
the condition that normals at a comer point remain normal during deformation, one 
obtains the desired stiffness properties . 

Since the heat shield is indeed a sandwich structure, it would appear that this 
element would furnish a superior representation to that given by the isotropic plate 
elements. In order to test this conjecture, analyses of a quadrant of the heat shield 
were performed for two gridworks as shown in Figure V-2 , a coarse mesh con- 
sisting of 28 sandwich elements and possessing 64 degrees of freedom, and a refined 
mesh consisting of 66 triangular sandwich elements and 156 degrees of freedom. 

This heat shield representation is assumed to have fixed support around the periphery 
along the bolt circle. 

In both analyses, the face thickness was 0.008 in. and the core thickness was 

0 

2.00 inches. The modulus of elasticity for the faces was 30 x 10 psi while the core 

4 

shear modulus, Gq, was taken as 1.8 x 10 psi. 

For comparison purposes, a refined gridwork analysis was performed utilizing 
quadrilateral and triangular plates, in the same manner as was described in SectionH. 
A comparison of the results of all three analyses is shown in Figure V-3, in the form 
of the centerline radial displacement profiles. 

Since the conventional plate representation is known to yield essentially con- 
vergent results for the heat shield for this degree of gridwork refinement, it is 
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A. Coarse Gridwork 



B. Fine Gridwork 

Figure V-2. Heat Shield Quadrant Triangular Sandwich Panel 
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apparent from Figure V-3 that a significant improvement is yet to be gained, through 
gridwork refinement, in the triangular sandwich element representation. On the other 
hand, the substantial improvement obtained by advancing the gridwork from coarse 
to fine indicates that the triangular sandwich element will be capable of yielding 
accurate results with sufficient refinement. 

The triangular sandwich element differs from the conventional plate element 
principally in the inclusion of a transverse shear flexibility. Consequently, the fore- 
going results prompted an examination of the relationship between core shear rigidity, 
G^, and gridwork size for a known solution. The problem chosen for examination was 

the simply- supported, uniformly loaded square plate shown in Figure V-4. The faces 

7 

are taken to be 0.005 in. with a modulus of elasticity of 10 psi. The core depth is 
assumed to be 1.0 in. 

Due to symmetry, only a quadrant of the plate need be analyzed. Three grid- 
works of elements have been utilized - 4, 16, and 64 elements. Analyses were per- 
formed for three values of core shear rigidity: 34,000 psi, 500,000 psi, and 3000 psi. 
The lowest value is the magnitude to be anticipated of conventional forms of sandwich 
construction. The high value is representative of an isotropic shape. 

The results of the analyses, in the form of the center deflection, are plotted as 
a function of gridwork refinement. Also shown are results obtained by means of the 
classical solutions derived in Reference 11. It is seen from these results that a good 
level of agreement exists for the weak core shear rigidity, while the agreement is 
very poor in the modeling of an isotropic plate. Agreement for the latter is improved 
by refinement of the gridwork but whether or not the refinement will lead to a solu- 
tion convergent upon the analytical solution is uncertain. Clearly, convergence will 
occur with only an inordinate degree of refinement. 

Thus, it is seen from the foregoing results that the triangular sandwich element 
furnishes workable results for a reasonable number of elements when the structure 
indeed evidences sandwich behavior. A sandwich form of construction may or may 
not evidence "sandwich behavior", i.e., significant shear deformation, depending upon 
the span, loadings, etc. 
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Deflection 
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The use of plate elements based on conventional isotropic plate flexure theory, 

i 

j which neglects shear deformation, has been justified by the view that the Apollo Aft 

Heat Shield has a relatively rigid core and a large span, so that the shear deforma- 
| tions are negligibly small when compared with those due to bending. In order to 

provide a quantitative estimate of the error produced by these assumptions, the case 

i 

| of a square sandwich plate was examined. The cross-sectional dimensions are 

taken to be similar to those of the Apollo Aft Heat Shield and are given in Figure V-5. 

I 

, Using curves generated in Reference 12 for rectangular sandwich panels, the maxi- 

mum (central) deflections for square panels of various sizes were obtained and com- 
| pared with the classical values for isotropic plates having the same bending stiffness. 

The percentage error is plotted in Figure V-5. 


For the heat shield, which is equivalent to a square plate of about 100 inches 
side length, the error is approximately 0.5%. 

__ / 

A convenient rule is to regard the error as negligible for L)p > 100 where 


\ = 


shear stiffness 


= flexural stiffness 


V) = panel width 

The curve of for various panel widths is plotted on Figure V-7. It can be 
seen that for ^ = 100 the error is approximately 2%. 

C. INTERELEMENT COMPATIBLE ELEMENTS 


As noted previously, rapid strides are being made in both the development of 
discrete element stiffness properties and in the definition of appropriate bases and 
concepts for the development of these properties. It is felt that these activities 
point towards idealization concepts which meet all of the needs of heat shield anal- 
ysis. The numerical evaluation and coding of these elements, however, has not 
advanced sufficiently far at this time to permit analysis of the Apollo Aft Heat Shield. 
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It is believed that discrete element idealizations for matrix displacement 
analysis should conform to the requirements of compatible displacement. If they 
meet these requirements, then it is assured that the displacement predictions will 
be smaller than the convergent solution and that the solution will monotronically 
approach the latter as the gridwork is refined in a consistent manner. Since discrete 
element formulations for matrix displacement analysis are generally based upon 
assumed displacement functions, there is no question regarding the satisfaction of 
compatibility within the elements ; the problem is one of so-contriving the respective 
element displacement functions so that the complete model satisfies compatibility. 
That is, interelement compatibility must be achieved. 


The element discussed previously - the triangular sandwich element - satisfies 
the interelement compatibility conditions. As shown, however, a very large number 
of these elements is necessary for an accurate representation. Thus, a more effi- 
cient approach is to utilize elements whose theoretical basis corresponds to that of 
conventional flexure theory, modified to account for the difference between mem- 
brane and flexural stiffness. Such elements are described in the following. 


From the standpoint of geometric representation, two geometric forms are 
essential - the triangle and the arbitrary quadrilateral. Advanced representations 
have recently (October 1965) been put forward for both of these forms. The advance- 
ment which the new quadrilateral and triangular plate flexure elements represent 
stems from a careful construction of the element displacement modes so as to 
satisfy completeness and geometric admissibility requirements with respect to the 
whole structure. Satisfaction of these requirements guarantees that the mathematical 
model represents a stiffer structure than the actual structure in the sense that it has 
a higher potential energy than the exact solution. Furthermore, as the number of 
degrees of freedom are increased (e.g., grid refinement) the solution will converge 
monotonically toward the exact solution. 

Suitable displacement behavior is described in the case of the quadrilateral 
element by defining a complete cubic polynomial mode shape over each of four zones 
(Figure V-6). 
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A total of ten coefficients are associated with each zone. Admissibility require- 
ments across zone interfaces reduce the total number of undetermined coefficients 
to sixteen. This is precisely the number needed to establish geometric admissibility 
across element interfaces. This is clear since the displacement along any line is 
cubic in a coordinate measured along the line and is, therefore, completely specified 
along the boundary, by the corner point displacements and slopes. Thus, displacement 
compatibility is established by the interelement matching of these quantities. The 
four undetermined coefficients which remain are needed to satisfy slope continuity 
across element interfaces. The slope along any line is quadratic in a coordinate 
measured along the line and requires, therefore, the specification of slope at one 
point along an element boundary in addition to the two comer points. The boundary 
midpoint are selected for computational convenience . The resulting quadrilateral 
flexure element has sixteen degrees of freedom manifest as one displacement and 
two rotations at each of the four corners and one rotation at the midpoint of each 
boundary. 

Completion of the formulation is effected in accordance with the well known 
Rayleigh-Ritz procedure. Substantial algebraic convenience is realized by the use 
of the oblique element coordinate shown in Figure V-6. Numerical evaluation of this 
advanced quadrilateral flexure element is being pursued at Bell Aerosystems 
Company. 

Complete, admissible displacement mode shapes for a triangular plate flexure 
element have been constructed using a procedure similar to that just described for 
the advanced quadrilateral element. As before, complete cubic polynomial mode 
shapes are assumed over zones of the element (Figure V-7). Continuity requirements 
across zone interfaces reduce the total number of undetermined coefficients to 
twelve. Arguments identical to those given for the quadrilateral element associate 
these coefficients with the nine corner point displacements and rotations and a 
rotation at each of the three boundary midpoints. This advanced triangular element 
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will, like the quadrilateral element, be incorporated into the Bell general structural 
analysis program. 

In conclusion, it is important to note that, not only are these triangular and 
quadrilateral flexure elements interelement compatible among other elements of 
their respective types, but quadrilateral elements may be interspersed with triangular 
elements without violation of interelement continuity requirements. 
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VI. ANALYTICAL SOLUTIONS 


A. GOVERNING EQUATIONS 

During the period when intensive efforts were directed toward the resolution of 
critical problems in the design of the Apollo heat shield for water impact conditions, 
the validity of the alternate paths toward numerical solutions for these problems was 
questioned. Thus, it becomes evident that classical solutions to simplified (but 
pertinent) problems were needed to verify that numerical analysis approaches. 
Furthermore, classical solutions can give a rapid, if approximate, measure of the 
quantities sought. The availability of convenient solutions will facilitate rapid analyses 
in the event new critical problems arise, and provide a means for preliminary design 
of future heat shields. With these considerations in mind, a portion of the work was 
devoted to the development of analytical solutions for a variety of heat shield pro- 
blems. For the Apollo heat shield water impact cases, the loading is not axisymmetric. 
Also, the effects of a sandwich type of structure must be included. As discussed 
elsewhere, the shear deformations of the sandwich core may be neglected (in view of 
the comparatively large span) and the sandwich solution is then similar to that for a 
uniform plate except that the membrane and bending stiffnesses differ and are repre- 
sented by the parameters and D* . 

With these considerations, the governing differential equations are 


vUtJ 


JL -2s. ( 

-*-) =o 

(VI-1) 
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where F is a stress function, , and the other symbols are defined 

Figure VI-1. These equations were derived by E. Reissner (References 


in 


The homogeneous solutions to Equations (VI-1) and (VI-2) are taken as: 


y ■ 0»i b«.r x Boa be>. % 


r A 03 + 

£ [B, ber x X -B a bei* X + A. X' ] <-ojm « 


T**l 


- - -IHIk 5 B oa ber X 1 - Bo i bev * + 

r\ “ D* ^ 


0, 


i [e il b«r„x *8, ^ a**"] 




(VI-3) 


(VI-4) 


It is seen that form < O there are basically four arbitrary constants to be deter- 
mined. These can be obtained from the following boundary conditions. 
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(a) For clamped edge 


no axial displacement 
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no radial slope 

c) u/ „ 

JV c “f 3s 

no circumferential strain 
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(b) For hinged edge 


(VI-5) 


(VI-6) 


(VI-7) 


(VI-8) 


Equations (VI-5), (VI-7), and (VI-8) remain while Equation (VI -6) is re- 
placed by the condition of no radial bending moment. This is in the form of Equation 
(VI-9). 


no radial bending moment 
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(VI-9) 


In the above equations the primes indicate that all the variables are expressed 
in nondimensional form. 


The development of the particular solutions for W and F has presented some 
problems which were not anticipated and which have necessitated the trial of a number 
of different approaches before a satisfactory solution could be obtained. 

Initially the loaded annular sector of Figure Vl-le was expressed in the form of 
a double Fourier Series, one of which was then transformed into a power series. 
Although a solution to the governing differential equations for this loading function 
could be obtained, comparisons of its degenerate form for the uniform load case with 
known solutions indicate that the solution was not valid. 


Eventually the difficulties were traced to the double Fourier, series which is 
apparently not valid on a shallow spherical shell. For a shallow shell Bessel Func- 
tions must be used and hence, the loading function has to be expressed through a 
Fourier-Bessel series. 
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This new solution is considerably more complex than the simple Fourier, series. 
Full details of the solution with results will be published as a supplement to this 
report. 

B. HEAT SHIELD WITH OVERHANG 


The analytical procedure presented above can be extended to obtain solutions to 
a heat shield with overhang. Such a case is defined in Figure (VI-3). The heat shield 
shown is supported along a circle of radius a. The portion of the shell inside the sup- 
porting circle is designated region "I" while the overhang portion is called region "H". 


The governing differential equations for this case remain the same as Equations 
(VI-1) and (VI-2). The homogeneous solutions to Equations (VI-1) and (VI-2) in Region 
(I) also remain of the same form as in Equations (VI-3) and (VI -4). In the overhang 
portion, however, no simplifications could be made except that one constant in 
may be disregarded due to the fact that the derivatives of ^ only are concerned. 
Therefore, the homogeneous solutions in Region (II) are 
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(VI-11) 

From Equations (VI-10) and (VI-11), it is seen that there are ten arbitrary 
constants for the case n = 0 and twelve constants for n 0 . 

To obtain the particular solutions, three loading cases are considered. These 
are (See Figure VI -2). 
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(1) The water impact area inside region (I) 

Yjl. o- 

Y<, — Ck 

(2) In both regions (I) and (II) 

Vi < 

To > a. 

(3) In the overhang portion only 

Yi > O 

Yi < r 0 ^ a* 

For the particular solution the same approach is used as for the shell without 
overhang. 

The boundary and continuity conditions needed for solutions are 
At r = a 

no axial displacements in both regions 


✓ 


UJ CD =o 

(VI-12) 

U^'n) r O 

(VI-13) 

circumferential strain in both regions 

( € » )x 

(VI-14) 


(VI-15) 

1 

meridional displacements in both regions 

Uu, --o 

(VI-16) 

U ci > “ ^ cn > 

(VI-17) 
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(VI-18) 


radial slopes are continuous 

t)r d v 

radial moments are continuous 

(.Hrt) x 1 (Vl-19) 


Atr=a h 

no radial bending moment 

(nr'r) u '-0 (VI-20) 

boundary shear must vanish 

(v/-£ (Vi-21) 

no radial membrane force 

(/Y r V) tt a ° (VI-22) 

no shear membrane force 

( n‘r 0 ) - O (VI-23) 

These above conditions are satisfied for the determination of the twelve 
arbitrary constants for the case n o. For n = o there are ten arbitrary constants, 
the conditions of Equations (VI-17) and (VI-23) are disregarded. 

The solution procedure from this point on is essentially the same as that 
described for the case without overhang. 

C. SOLUTION FOR UNIFORM THERMAL FORCE AND UNIFORM THERMAL 
MOMENT 

Under the conditions of uniform thermal loading the governing differential 
equations (Equations VI-1) and (VI-2) reduce to the homogeneous form 

^2- V*w'-0 (VI-24) 

v r D 5 

- STl F ’ * ° (VI-25) 

Since the current problem is axisymmetric the homogeneous solutions can be written 
as 
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(VI-26) 


W = C, be* K - C a belt XV tj 
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(VI-27) 


It is to be noted that, as yet, the effects of the thermal loading have not been 
introduced. In the case of uniform thermal loading, stresses are only introduced by 
restraining the thermal deformation and hence the magnitude of the thermal loading 
appear only when the boundary conditions are considered. The shell is fully fixed at 
the bolt circle. The boundary conditions are then 


W - O 
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where the circumferential strain is given by 
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and /V . is the thermal force. 


It is to be noted, as was to have been anticipated, that the thermal moment 


tern* does not appear explicitly in this case. 


The three relevant constants c^ C2, c 3 are then obtained by the solution of the 
simultaneous equations as 
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For the case under consideration the following numerical values were used. 

<K $- lU.t , O s - 3 X IO fc , Pw» s O. 4 4 , Ok* 5*7.0. 

h/V,* s 1 9 of I V>. , No*. s 3 H 2 (» lb. / ik. 

Th**s i- ~ *4* - 14. 98 3 , = O, 3 
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I 

I 

i 

c, Z ».4.* * O A IO~ 3 , C- a. - 2. 3SS A io' A J Cj - 7. 1C Ti A to A 


Central Deflection = J?(c,+-c 3 )= | . 0 8 8 =2 in. 

A complete profile of the displacements along the axis of symmetry is given in 
Figure 

At center of heat shield 

(-7- beiVX^ - 57.«*G IkV/n. 

M V - - -2i c* • We = - 1 4T» + 1 ro*»- * 2.C* iw /* *•- lit i. 


For the chosen section, the details are shown in the sketch. 




The position of the neutral axis is 
1.953 from the lower surface. 


0.008 in. 


T 

I 2.0 in. 


:0.008 in. 


Using the calculated membrane force 
H t = 57.66 lb/in., the stresses in the 
ablator and steel faces are 

(TXi fcfc - I 3.7° P* 4 (Compressive) 
(Tia* = H O 4 . \ f * i (Compressive) 


The bending moment M *• = 26 lb/ in./ in. produces the following stresses . 


In lower steel face (J* s - \ O A SM (Tensile) 

In upper edge of ablator (T$ - S A . O ps* (Compressive) 
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VH. CONCLUDING REMARKS 


The discrete element approach to structural analysis has been shown to be an 
appropriate means for the analysis of the Apollo aft heat shield where the seemingly 
regular geometry would appear to make attractive classical solution approaches. In 
the actual circumstances, the analysis must account for irregular loadings, tempera- 
tures, and a wide variety of possible support conditions. More significant perhaps is 
the irregularity of the cross-sectional dimensions as produced by scalloping of the 
sandwich faces and the variable thickness ablator, and by the portion of the heat shield 
which overhangs the ablator. All of these effects are treated directly by means of 
discrete elements of triangular and quadrilateral form. 

Further developmental work in connection with discrete element properties is 
nevertheless indicated. Two promising directions have been found. In one, sandwich- 
type elements can be employed. Presently, this approach required implementation with 
respect to the refinement of the idealization, wherein a very much larger number of 
degrees of freedom must be utilized. Also, the available element relationships which 
are derived on the basis of equal face thicknesses, must be extended to account for a 
load carrying ablator attached at one face. 

It should be additionally noted that if design requirements involve the determina- 
tion of peaked stresses, an approach based on finite elements (or finite differences) 
will be inadequate, or in the least inefficient, if the problem is approached through use 
of a single idealization of the complete heat shield. Logically, a gridwork of reasonable 
refinement should be employed to obtain the overall stress and displacement pattern. 
Then, a high degree of refinement should be employed in the idealization of an isolated 
portion of the structure, where the determination of peaked stresses is desired. Pro- 
cedures for this type of approach are easily defined, but have not been properly tested 
in practice. 

With respect to buckling, a finite deformation -finite displacement formulation 
would be a logical step toward improved prediction of behavior of the heat shield. The 
present instability analysis capability exhibits complete membrane -bending coupling 
between discrete elements of an assembly. With reference to a single discrete element, 
however, only the effect of membrane action on bending stiffness is considered. This 
latter procedure is a conceptual extension of the linear formulation which allows the 
prediction of critical loads for structures which experience small deformations and 
small displacements prior to buckling. 

Although the finite displacement formulation involves major modification of the 
instability analysis capability, it could be expected to significantly improve the predic- 
tion of behavior. As stated by B. Budiansky and J. Sanders as a result of their 
examination of certain phases of this Apollo heat shield structural analyses (Reference 
18) "should a static nonlinear analysis be done? On the basis of some estimates we 
made we concluded that one should be done." In this same report, emphasis is given 
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to the view that since a nonlinear analysis is feasible,it is important to perform it and 
dispense with the empirically -modified linear result, since the empirical factor is of 
ill -defined reliability. 

As a final comment, the determination of thermal stresses in the heat shield is 
a subject requiring more extensive and detailed attention that it could be given in this 
effort. Classical solutions for significant thermal stress situations in shallow shells 
are virtually nonexistent. This being the case, proper verification of the adequacy and 
special safeguards to be taken in a discrete element approach is difficult to achieve. 

In a specific, related, problem the extremely high predicted thermal stresses 
for the heat shield are indicative of an unrealistic definition of the support conditions. 
The elastic nature of the support in the radial direction in the plane of the bolt circle 
should be taken into account. Also, no recognition has been given to the existence 
of initial assembly stresses. It is known that the ablator is applied to the heat shield 
at a temperature other than ambient room temperature. Thus, at room temperature, 
thermal stresses are produced in consequence of the differences in expansion coef- 
ficient of the ablator and the sandwich. 
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APPENDIX A 


MATRIX DISCRETE ELEMENT ANALYSIS PROCEDURE - 
LINEAR STATIC ANALYSIS 


1. GENERAL CONCEPTS 

The purpose of this appendix is to develop, from basic principles, the concepts 
of elastic linear, discrete element structural analysis (10)(16), as applied to static 
loading conditions. 


In this method of analysis, structures are idealized as systems of connected 
discrete elements. The points of connection are called "reference" or "node" points. 
Each class of discrete element (bar, triangular plate, beam segment, etc.) possesses 
a finite number of connection points, the specific number in a given case being 
dictated by the number of parameters needed to define the variation of the edge stress 
systems acting upon the element. A hypothetical element, a rectangular plate seg- 
ment of a plate component, having four reference points, is shown below. (Figure A-l) 



Y 



For any such element, it is first necessary to derive relationships between the 
displacements {& e \ of the boundary points and the forces l F \ acting at these 
points. The node point forces are statically equivalent to the stresses that actually 
exist on the edge areas subtended by the point. On the basis of assumptions as to 
element deformational behavior, it is possible to establish the desired relationships 
in matrix form, as 

M- : [K]KWF'j 
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where 

[KJ is the "element stiffness matrix”, 

(M are the displacements of the points on the element (node points) which join 
with the adjacent elements or with the supports, 

| P ^ are the stress resultants, or equivalent forces, at the element node points, 

l F m j are the "thermal forces", representing the effect of a temperature change 
with the element. 

Methods for the formulation of the terms of this equation are detailed in 
References 11 and 16. 


Once the element relationships have been evaluated, the elements are assembled 
to form the complete analytical model of the structure by joining all elements at their 
respective juncture points and applying, in the process, the requirements of juncture 
point equilibrium and compatibility. The theoretical basis for this operation is as 
follows. 

The components of internal loads ^ F $ and net external loads { P ^ at each point 
are related by equilibrium requirements; i.e., ft, etc. The respective coordinate 
displacements of the corner points of all elements meeting at a point are equal, a 
requirement that satisfies compatibility. The result is that the stiffness matrix [Klfor 
the complete structure can be assembled by merely adding element stiffness coefficients 
having identical subscripts. This yields a set of algebraic equations: 


H- Ck] ^ a] + [P“J 


(A-2) 


The matrix L K 3 will henceforth be referred to as the "master" stiffness matrix. 
Displacement boundary conditions can be readily imposed by assigning the pertinent 
A 's their known values (usually zero). The matrix [Kjwill be altered in the process, 
and, taking note of this by utilizing the subscript R, the solution to the altered 
Equation (A -3) becomes (if matrix inversion is utilized), 




(A -3) 


where is] represents the set of displacement influence coefficients. This is the 
equilibrium static solution. The displacement influence coefficients are directly use- 
ful for subsequent dynamic analyses. 

To obtain the stresses from the displacement solution, the displacement vectors 
for the respective elements [ are first selected from the total column of displace- 
ments. Then, each such vector is multiplied by the associated stiffness matrix (Equa- 
tion A-l) to determine the node point forces j F* e J . In an additional step, the node 
point forces can be transformed into the corresponding stresses. It has been found 
more efficient, however, to form at the outset direct relationships between the element 
stresses and the node point displacements, as follows: 
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*l<r~i (A-4) 

where ] are the stress values which characterize the distribution of stress within 
an element, fa^jare the stresses corresponding to the thermal forces, and tS lis 
known as the "element stress matrix". The procedure followed, therefore, is to 
establish, first, the stress matrices at the start of a computation. When the displace- 
ment vectors for the respective elements | A e 3 are evaluated, they are premultiplied 
by the corresponding element stress matrices to obtain the solutions for stress. 

2. BELL GENERAL PURPOSE STRUCTURAL ANALYSIS COMPUTER PROGRAM 

The Bell "General Purpose" Matrix Structural Analysis Program, coded for 
operation on an IBM 7090 computer, is designed to accept the basic information des- 
cribing a problem, establish its formulation, and perform all computations required 
for the development of the desired stress and/or displacement results. The opera- 
tional procedures for the program conform to the theory outlined earlier. 

The basic problem information (input to the program) consists of: 

(1) Dimensions of the structure 

(2) Load and temperature conditions of interest 

(3) Material mechanical properties 

(4) Operational controls; i.e., specification of the desired printout items, etc. 

(5) Designation of the discrete elements 

From an operational standpoint, the Bell General Purpose Structural Analysis 
Program consists of three major computational routines: 

(1) A library of element stiffness relationships 

(2) A routine wherein the master stiffness matrix is calculated, boundary 
conditions are applied, and the matrix inverted. 

(3) A routine which selects information from (1) and the resulting inverse 
from (2), and calculates stress, displacement, etc. 

Routine (1) is the key to the versatility of the program, since the capability to 
analyze a given type of configuration is dependent upon the availability, in the element 
library, of relationships for elements of the proper geometric form and behavior. A 
listing of the more commonly used elements in the program is given as Table A-l. 
These elements are pictured in Figure A -2. Elements employed in the subject analyses 
are discussed in the next section. 
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Three of the detailed operational capabilities of the program are of specific 
interest; 

(1) Although element properties are stored in the program in algebraic form 
and are automatically evaluated on the basis of input data, it is also per- 
missible to input directly an evaluated element stiffness matrix. Thus, if 
the stiffness coefficients for a major component have been evaluated 
elsewhere, this component can be represented in the analysis as a single 
discrete element. 

(2) Special coordinate axes can be established at any or all points. The force 
and displacement vectors at each reference point are initially referenced 
to a single set of axes for the complete structure (system axes) . The need 
for utilization of this capability arises, for example, in the following cases: 

(a) A structure is constrained to displace in directions other than the sys- 
tem axes (Figure A-3a). Thus, thex'-y' axes must be defined at all the 
support points and used to specify the y' -direction constraints. 

(b) On a shell in flexure, when the individual elements are each flat plates, 
the two-dimensional moment vectors on the respective elements 
meeting at a point produce on extremely small net vector in the direc- 
tion normal to the tangent plane at the point (Figure A -3b). If the 
coordinate axes are made to correspond to the tangent plane , this 
degree -of-freedom normal to the tangent plane can be suppressed with 
a consequent improvement in the accuracy of the solution. 

(3) The inverted stiffness matrix (i.e., the flexibility matrix) can be stored on 
' tape and any or all elements of this matrix can be recalled at a later time 

for the purpose of subsequent dynamic analyses. 
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APPENDIX B 

DYNAMIC RESPONSE ANALYSIS PROCEDURE 


1. GENERAL 

The equations of motion for a system of grid points , describing an elastic struc- 
ture, can be written in the following matrix form: 

[M3 l&l + CD3 [Aj + flO/A] = l (B-l) 

where 


CM3 
L 03 
[K3 


is the mass matrix 

is the damping matrix 

is the pertinent stiffness matrix 

are displacement, velocity and acceleration vectors, respectively 
is a forcing function matrix 


The displacement vector £ A J can represent motion along and about coordinate axes. 
The mass, damping and stiffness matrices will then consist of terms consistent with 
these motions. 


Equation (B-l) is applicable to any elastic structure. The component matrices, 
however, depend on the coordinate system employed. The displacement degrees of 
freedom can be identical to those of the structural idealization grid system or 

can be derived from the natural modes of vibration of the structure. The coordinates 
resulting from the latter approach are called generalized coordinates and are dis- 
cussed in detail in this Appendix. The appendix is concluded with a description of the 
integration of Equation (B-l). 

The use of generalized coordinates (or the so-called normal mode approach) re- 
sults in equations of motion with the following characteristics: 

(a) The first eight to ten modes of vibration will generally be sufficient for 
dynamic response analysis purposes, as compared to perhaps 300-500 
degrees of freedom required for the accurate determination of the flexibility 
characteristics when the structural grid system is used. Therefore, lower 
order matrices in Equation (B-l) will result. 

(b) The mass matrix C M J will be a diagonal matrix as there will be no mass 
coupling between the natural modes . 

(c) The stiffness matrix [ K 3 will be a diagonal matrix since there will be no 
elastic coupling between natural modes. 
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The adoption of the model approach in the dynamic response analysis portion of 
the investigation permits the initial determination of stiffness characteristics on the 
refined basis described previously, while enabling the performance of computationally 
efficient analyses. 

2. DYNAMIC RESPONSE EQUATION USING GENERALIZED COORDINATES 

Consider the structural grid system sketched below where the ith grid point 
has an assigned mass, mi, and is undergoing a displacement wj. In this instance, 
the displacement vector 3 consists of the single translational degree of freedom 
w at each grid point. The following results can readily be extended to include other 
degrees of freedom. 



Figure B-l. Typical Grid System 


It is assumed that the displacement at point i, due to a dynamic load F( x- ^ tj* } t ) 
can be represented by the series 

'UTj ft) % j(t) (B_2) 

where /wj ( j <j;hs the ^ ^ normalized vibration mode shape evaluated at the i ~ 
point and (t) is the generalized coordinate in the j vibration mode. 
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t 


or 


*1A 

The kinetic energy at the L point may be written as 

Ti * 5 [»y ( *•, ft, * ) J 

I- * k <*.■ C| ( ">• ?,•]* 

Expansion of Equation (B-4) yields 

X = i 2cr*i, Z +k sm; iz 


+ - • • 


+ i i. it + i *7 &3 i, i 3 + * * * 

Since natural modes are orthogonal the sum of terms of the form 

I • « / | 

2 'm; j I *'>* ^ j 7 = * 

Thus, Equation (B-5) becomes 

Ti = £ '»»** %z + *5* 4/ + * * * 4 ^ J 

for j = k. 


(B-3) 


(B-4) 


(B-5) 


(B-6) 


j m t \ 

The kinetic energy term in Lagrange's equation of motion is given by ^ 1 

where T is the total kinetic energy of the system, i.e. $ f • 


T= I TV = Z 2. /77ij [ ~urf- £+ + - •*+ ^ + - - ] 


(B-7) 


Thus, we can write 

c£ 


J * 


) ) j 


etc. 


«£ / £E ) r £ ^fi 
£, / *»/ 


(B-8) 


or more generally 



where Mi s £ /wy ft)J is the generalized mass. 

* isl 


(B-9) 
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The strain energy of the system is derivable by analogy to a simple spring- 
mass system vibrating in a natural mode. When the entire grid system is vibrating in 
a natural mode it can be said to have a mass Mf (the generalized mass), displace- 
ment § j.' (the generalized displacement), stiffness K j (the generalized stiffness) and 
natural frequency uj^ . The net potential energy stored by this system is given by 
(see Figure B-2). 


U.. 


^Z ZZZ /Z ZZZZ. 


'8 st 


1 

| £ j 

M. 

J 

g 




Figure B-2.. Displacement 


(B-10) 


Ui - {£K a -(S ir +S 4 -AiK/S4]-(f )€,* 

r ^ j 

Since Equation (B-10) can be rewritten as 

V* r ? i 

Now wj- Hj 

hence . I u>f . 

The total potential energy of the grid system is given by the following expression: 

U = 2 Uj - 2 i U>1 + *" + i <“i (B ' 14) 

V 1 5 4 


(B-ll) 

(B-12) 

(B-13) 


Equation (B-14) is a quadratic function of the generalized coordinates and agrees 
with the approach discussed in Reference 
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The potential energy term in Lagrange’s equation of motion is given by 
where V is the total potential energy of the system. Thus ° i 

^Si (B- 

g - H,* | z 


or more generally 


c)T7 . 


« Kj €. 


where K j, = is the generalized stiffness and OJ^' is the frequency of the 

vibration mode . 

In accordance with Reference 17 an energy dissipation function associated with 
structural damping can be setup as follows 


7s . _ J_ u)$ 

* Z u> 


H . 


is the damping term associated with the mode, fj. is the damping coefficient 
in the f ~ mode and to is the response frequency. The damping term in Lagrange's 

equations of motion is given bv where J> is the total energy dissipated. Thus, 
we may write 

S D . Dj i; (B-17) 

where Dj,* MjOf the generalized damping. 
co 

The total virtual work done on the system due to an incremental displacement 
(or virtual displacement) in each degree of freedom is 


- 21 Su- ( *<: a qi 3 -t ) 


Substitution of Equation (B-2) into Equation (B-18) yields 


OD V 


5 W e - ^ ^ * c •> ~* L ' o ^ ) ) & £ £ 
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The generalized force term in Lagrange's equations of motion is given by 

HiVk) 


<2 


thus 





(B-21) 


Equations (B-9), (B-16), (B-17), and (B-21) are summed to form the differential 
equations of motion of the system. 


or in matrix form 


(B-22) 


(B-23) 


Equation (B-23) is identical in form to Equation (B-l) with the exception that 
the matrix coefficients are diagonal; that is, the equations of motion are uncoupled. 

For a finite number of vibration modes, say j = 3, a finite number of grid 
points, say N = 6, and forcing functions acting only at grid points 1 and 2, three 
equations in three unknowns are obtained as follows: 
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Note that the number of modes determines the matrix order of the left-hand 
side of the equation while the number of modes and number of forcing functions de- 
termine the matrix order on the right-hand side. 

It is a simple matter to determine the true displacement at a grid point by use 
of Equation (B-2) once the generalized displacements as a function of time are known. 
Suppose that in the present illustrative example displacements at grid points 1, 3, 5, 
and 6 are desired. Use of Equation (B-2) yields the matrix expression. 



1 

X 

■ft 

K 

1 




* r t(x s ,‘l f ) 





(B-25) 


This same matrix expression is used to determine grid point velocities and accelera- 
tions by replacing the CO j matrix by^£ (O j and | ^ Ct) J matrices. 

3. INTEGRATION OF THE EQUATIONS OF MOTION 

The evaluation of dynamic response, in the form of time histories for A , A 
, * • 

and A involves the determination of instantaneous acceleration over a small interval 
of time , followed by the performance of single and double integrations to obtain velocity 
and displacement changes over the same interval. The velocity and displacement 
changes are then added to the initial starting values , and the process is repeated for 
the next interval. This process forms the basis of the Bell-Matrix-Iteration-Time- 
History program which is described below. Note that the method applies equally well 
to Equation (B-l) or (B-23). 

The method proceeds as follows. First, Equation (B-l) is written in the form 
of a solution for m 


/a]- CM j" {{Fj -[0]JAj -[KjfAjJ 


(B-26) 


The matrices^FjjfMjjliD JmdC Klare known. For any starting time t G , the values of 
^Ajand ^A j are prescribed. The latter are zero if the system is initially at rest. 

The foregoing values of acceleration computed for t Q are assumed to exist over 
a small time interval, t^ - t Q , and values of ^A | are computed from 

MJ.* i*'] (*•'*•)({« ],+ /*}.} (B-27) 
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which will be recognized as a simple time integration. Then, the displacement at the 
end of ti is then computed from a second integration 

(B-28) 


It will be noted that average velocities and displacements are obtained simply as one- 
half the sum of the initial and final values . While much more elaborate formulations 
could be used, they are considered less desirable than using closely spaced time 
intervals . 

It will be noted that some error will exist in the above computations when the 
system is elastically deformed since the deflection pattern at t^ will not be the same 
as at t D , causing accelerations at t Q to be different than accelerations at t^ via 
Equation (B-26). To eliminate this error, an iterative procedure is used. This in- 
volves recomputing accelerations for t^ using the velocities and displacements from 
the initial computation. These accelerations are then used in Equation (B-27) to ob- 
tain more accurate average velocities than were obtained the first time when it was 
assumed that the initial acceleration remained constant over the time interval. The 
new velocities A, are then used in Equation (B-28) to obtain more accurate final 
displacements . 

The cycle of acceleration, velocity and displacement computations is repeated 
until convergence to a satisfactory number of digits. Convergence is extremely 
rapid, particularly when time intervals are small enough to avoid rapid changes in 
acceleration between time intervals. After convergence has been obtained, the com- 
putation proceeds to t = t 2 and repeats the process already described here for t = t±, 
etc. In this manner, time histories for A , A and A are computed. 
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APPENDIX C 


MATRIX INSTABILITY ANALYSIS PROCEDURE- LINEAR 


1. THEORETICAL BASIS 

Elastic instability can occur for prismatic thin-walled structures or thin 
shells, where the behavior across the thickness can be subdivided into "flexural" 
and "midplane" behavior. By virtue of displacements normal to the midplane, the 
membrane forces have components which tend to enhance these displacements. 

When their magnitude is sufficiently large they produce — in the context of linear 
theory - infinitely large displacements for infinitesimal changes of the applied load. 

The values of applied load which produce the membrane stress state to cause this 
instability are the "critical loads". 

In the case of straight beams and flat plates, the membrane forces are entirely 
independent of flexural behavior if finite displacements are disregarded in the formu- 
lation of the strain-displacement equations. For more general configurations such 
as shells, however, the applied loads produce membrane and flexural deformations 
which are related to each other. These influences are accounted for in the present 
elastic instability analysis process , although the nonlinear terms in the strain-dis- 
placement equations continue to be disregarded. 

2. THEORETICAL BASIS 

Although the interaction between the membrane and flexural behaviors is 
accounted for in the overall analysis procedure for the assembled structure, the 
elements are in themselves flat elements and their membrane and flexural stiffness 
relationships can be developed separately. 

In view of this independence of element flexural and membrane behaviors, one 
can partition a plate element (Figure C-l) stiffness matrix, given earlier as Equation 
A-l of Appendix A, as follows 
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where the primes designate "local coordinates (the coordinate system affixed to the 
element) . 
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The membrane stiffness for the flat plate element is not directly affected by 
instability considerations. As shown subsequently, however, element stiffness for 
flexural behavior in the presence of membrane forces becomes the sum of two com- 
ponent stiffnesses: 

[K,]-[K b ] ♦[*] < c - 2 > 

where C^blis the stiffness for flexural behavior in the absence of instability effects 
(the conventional flexural stiffness matrix) and [ft ] represents the effects of the 
element membrane forces on the element flexural behavior. The terms of CnJ con- 
sist of the dimensions of the element and the values of the membrane forces acting 
upon the element. Material properties do not appear in the fnj matrix. 


Using Equation (C-2), the total element stiffness matrix can be written as 
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(C-3) 


The coupling of membrane and flexural behaviors in a complete analysis stems 
from the coupling of the elements at the joints. In the general case, the respective 
elements will be randomly oriented with respect to the system axes. When joining 
the elements at the node points (i.e. , when applying the conditions of joint equili- 
brium) , a given force component in a particular lpcal direction will be transformed 
into components in all three system directions. As shown in Figure C-2 in two 
dimensions, the Fx A force is transformed into components of both Fx A andF / A 
when the equations of equilibrium are written at point i. By virtue of this operation, 
the respective force components of the elements in local coordinate directions con- 
tribute to all three directions in the system-coordinate equations for the complete 
structure. Thus, in the presence of instability effects, the latter have the following 
general form: 
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Equations (C-4) are the governing equations for the total problem. In parti- 
cular, it is to be noted that the evaluation of the [.Nj matrix requires a knowledge of 
the membrane stress state. Classical analysis techniques commonly assume that 
the membrane stress state is unaffected by flexural deformation. Based on this 
assumption, a membrane stress state is readily determined. The same approach 
can, of course, be taken on a discrete element analysis. The determinate membrane 
stress state would be hand computed and used to calculate the [N J matrix. 

The digital computation procedure has been arranged, however, so that the 
latter simplification need not be made. Consider first an equilibrium analysis, where 
the applied loads are not large enough to cause buckling. The analysis would begin 
with the^Nj matrix first set equal to zero - instability effects are disregarded. The 
solution for displacements is then 

M3 = [K] ' JP, M } (c-5) 

The so-determined values of Ml are next used in the element stress equations (see 
Equation A-4 of Appendix A and the related discussion) to obtain a first guess for 
the membrane stresses. These membrane stresses account for the influence of 
bending. With these membrane stress values, the matrix [ N/ J can be constructed 
and Equation (C-4) is solved as follows 

\b]- f[KMN]]''[P M] < c - 6 > 

Again, element membrane stresses are evaluated and the process continues iteratively 
until convergence. 

To predict buckling, the approach taken involves a sequence of eigenvalue 
determinations for a succession of assumed applied load levels. In a given eigen- 
value determination, based on an assumed applied load level, it is assumed that all 
midplane forces are at a fixed ratio to one another at all levels of applied load — 
from the load level chosen up to the specific value predicted to cause instability in 
the eigenvalue determination. Thus, Equation (C-4) is written as: 

= [ K ] fa j + A [n][aJ (C-7) 

where [N] is constructed from membrane stresses consistent with the assumed 
applied load level and A is the scalar multiplier which defines the buckling level con- 
sistent with the membrane stress state for the assumed applied load level. Setting 
\p t M } equal to zero for the purposes of eigenvalue analysis 

o * + A O ]{a] 

X hi - -OJ"["]H 

Using matrix iteration, the above can be solved for the eigenvalues A* and the 
associated eigenvectors | A;]. There will be as many such eigenvalues as there are 


(C-8) 

(C-9) 
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I 

equations in (C-8), but the only eigenvalue of interest is the largest value of y re- 
presenting the smallest \ and, therefore, the lowest magnitude load at which elastic 
instability will be experienced. 

The actual buckling load will, of course, occur only when the determined A = 
1.0. To accomplish this, atrial and error process is attempted. First, an obviously 
low value of applied load is assumed, resulting in a high value of ( > 1.0). This is 
shown as Point A in Figure C-3. Correspondingly, if the assumed load level is higher 
than that which causes buckling, the predicted A will be less than 1.0 (Point B). By 
obtaining values of A at various assumed load levels (Points C & D) , the load level 
for buckling can be determined (Point E) . 

It should be noted that the determination of A for an assumed buckling load 
level is accomplished only after a sequence of iterative analyses has been performed 
(in the manner described previously) to determine a convergent membrane stress 
system and, therefore, INI matrix. 

3. ELEMENT RELATIONSHIPS 


Procedures for formulating discrete element force-displacement equations 
were reviewed earlier, in Appendix A. These procedures are based on the use of 
Castigliano's Theorem, which, in the formulation of a stiffness matrix takes the 
form „ 
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(C-10) 


where LB] is the matrix defining the element node point displacements \ j as a 
function of the undetermined parameters, of the assumed displacement functions for 
this element, i.e. 

M < c - n > 

and [c 7 is a matrix, each row of which represents the derivative of the strain 
energy of deformation of the element (U) with respect to one of the undetermined 
parameters. 


Reference 6 shows that the above approach can be employed in the derivation 
of stiffness matrices which include instability effects if [C lis expressed as the sum 
of two matrices 


[C] -- [ Cj+[CJ (C-12) 

where F ] derives from the strain energy for flexure along, ( Ut, ) and l C-nl 
represents the work done by the edge membrane forces acting through the edge dis- 
placements of the element (U n ). 
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Typical Node Point 


Figure C-l. Flat Plate Element 




(a) Local Coordinate Forces 
at Point i 

(Direct Forces Only) 


(b) System Coordinate 
Forces at Point i 


Figure C-2. Joint Equilibrium Conditions 



Figure C-3. Iterative Determination of Critical Load 
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For plates, 



l/ b was given previously, in Appendix A 
Substituting (C-12) and (C-10) 

[k]' (LbT‘) [C b lt Bj'+ ( t&i' ? i c IB J' <c-!4, 

so that, by comparison with Equation (C-2) 

[Ki ■ { 1 6J") T r v [b]” 

In] . ( l'bjT [ c J [ B 

Using the above concepts, and the assumed plate element deflectional nodes 
described in Appendix A, detailed derivations of the flexural stiffnesses for the 
triangle and quadrilateral have been performed and are presented in Reference 3. 
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